arXiv: 1507.00012v2 [astro-ph.GA] 26 Feb 2016 


MNRAS 000, 1-?? (2015) 


Preprint 1 March 2016 


Compiled using MNRAS IAI^X style file v3.0 


Simulator of GAlaxy Millimetre/submillimetre Emission 
(siGAME*) : CO emission from massive z — 2 main-sequence 
galaxies 

Karen P. Olsen , 1 f, Thomas R. Greve, 2 f Christian Brinch, 3,4 f Jesper Sommer-Larsen, 
Jesper Rasmussen , 1,7 Sune Toft , 1 Andrew Zirm 1 

1 Dark Cosmology Centre, Niels Bohr Institute, University of Copenhagen, Juliane Maries Vej 30, DK-2100 Copenhagen, Denmark 
2 Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK 
3 Centre for Star and Planet Formation (Starplan) and Niels Bohr Institute, University of Copenhagen, 

Juliane Maries Vej 30, DK-2100 Copenhagen, Denmark 

4 DeIC, Technical University of Denmark, Building 309, DK-2800 Kgs. Lyngby, Denmark 
5 Excellence Cluster Universe, Boltzmannstr. 2, 85748 Garching, Germany 
6 Marie Kruses Skole, Stavnsholtvej 29-31, DK-3520 Farum, Denmark 

' Department of Physics, Technical University of Denmark, Building 309, DK-2800 Kgs. Lyngby, Denmark 


ABSTRACT 

We present SIGAME (Simulator of GAlaxy Millimetre/submillimetre Emission), a new 
numerical code designed to simulate the 12 CO rotational line spectrum of galaxies. 
Using sub-grid physics recipes to post-process the outputs of smoothed particle hy¬ 
drodynamics (SPH) simulations, a molecular gas phase is condensed out of the hot 
and partly ionised SPH gas. The gas is subjected to far-UV radiation fields and CR 
ionisation rates which are set to scale with the local star formation rate volume density. 
Level populations and radiative transport of the CO lines are solved with the 3-D ra¬ 
diative transfer code LIME. We have applied SIGAME to cosmological SPH simulations 
of three disk galaxies at z = 2 with stellar masses in the range ~ 0.5 — 2xlO n M 0 and 
star formation rates ~ 40 — 140 M 0 yr -1 . Global CO luminosities and line ratios are 
in agreement with observations of disk galaxies at z ~ 2 up to and including J = 3 — 2 
but falling short of the few existing J = 5 — 4 observations. The central 5 kpc regions 
of our galaxies have CO 3 — 2/1 — 0 and 7 — 6/1 — 0 brightness temperature ratios 
of ~ 0.55 — 0.65 and ~ 0.02 — 0.08, respectively, while further out in the disk the 
ratios drop to more quiescent values of ~ 0.5 and < 0.01. Global CO-to-H 2 conver¬ 
sion (aco) factors are ~ 1.5M 0 pc -2 (Kkms -1 ) -1 , i.e., ~ 2 — 3x below the typically 
adopted values for disk galaxies, and aco increases with radius, in agreement with 
observations of nearby galaxies. Adopting a top-heavy Giant Molecular Cloud (GMC) 
mass spectrum does not significantly change the results. Steepening the GMC density 
profiles leads to higher global line ratios for J up > 3 and CO-to-H 2 conversion factors 
(~ 3.6M 0 pc -2 (Kkms -1 ) -1 ). 

Key words: radiative transfer - methods: numerical - ISM: clouds - ISM: lines and 
bands - galaxies: high-redshift - galaxies: ISM 


1 INTRODUCTION 

Inferring the physical properties of the molecular gas in 
galaxies is an important prerequisite for understanding their 
star formation. Carbon monoxide (CO) has proven to be a 
valuable tracer of the molecular interstellar medium (ISM), 

* SIGAME means ‘follow me’ in Spanish. 

f Email: karen@dark-cosmology.dk (KPO); t.greve@ucl.ac.uk 
(TRG); brinch@nbi.dk (CB) 


its rotational transitions being sensitive to a range of gas 
densities and temperatures through collisional excitation by 
Ht. Observations of CO lines have therefore been used ex¬ 
tensively to probe the conditions of the molecular ISM in 
galaxies at high and low redshifts. However, the number of 
z > 1 galaxies detected in CO is still relatively modest - 
little more than about 200 at the time of writing - with 
only a fraction of those having multiple CO transitions ob¬ 
served (see review by Carilli & Walter 2013). Furthermore, 
the vast majority of high-z CO measurements are of extreme 
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objects such as submillimetre-selected galaxies (SMGs, see 
Frayer et al. 1998; Neri et al. 2003; Bothwell et al. 2013; 
Yun et al. 2015) and powerful QSOs with infra-red luminosi¬ 
ties in excess of 10 12 L@ (Barvainis et al. 1994; Weifi et al. 
2007; Riechers et al. 2011). Very few normal star forming 
galaxies at high redshifts, i.e., galaxies on the galaxy ‘main 
sequence’ (MS), have been observed in CO to date, and as 
a result little is known, however, about the ISM conditions 
of normal star-forming galaxies at high redshifts. Some of 
the best studied cases are the four BzK-selected galaxies 
BzK—4171, BzK—16000, and BzK—21000 at 2 ~ 1.5, where 
up to four CO transitions, including CO (1—0), have been 
observed (Dannerbauer et al. 2009; Aravena et al. 2010; 
Daddi et al. 2015). While the initial analysis of the low- 
J {Jy -ip < 3) transitions in these galaxies suggested a quies¬ 
cent ISM reminiscent of the Milky Way (MW) (Dannerbauer 
et al. 2009), the more recent CO(5 —4) observations require 
a second much denser, and likely warmer, ISM component, 
which is thought to be directly associated with star forming 
clumps within the galaxies (Daddi et al. 2015). Studies such 
as these illustrate the importance of sampling both the low- 
and high-J parts of the CO Spectral Line Energy Distribu¬ 
tion (SLED) in order to obtain a comprehensive picture of 
the molecular gas in galaxies. 

With the Atacama Large Millimetre/submillimetre Ar¬ 
ray (ALMA) now in operation, the collection of CO data is 
gaining rapid pace as is our understanding of the molecular 
gas in distant galaxies. Galaxy evolution simulations that in¬ 
clude molecular line emission (such as CO) are an important 
tool with which to interpret the observed molecular emission 
of galaxies caught in various evolutionary stages. The past 
decade saw the development of the first detailed simulations 
of CO line emission from galaxies, which consisted of SPH 
galaxy simulations (Narayanan et al. 2006, 2011; Narayanan 
& Krumholz 2014; Greve & Sommer-Larsen 2008), or semi- 
analytical models (Obreschkow et al. 2009, 2011; Lagos 
et al. 2012; Munoz & Furlanetto 2013; Popping et al. 2014), 
combined with sub-grid physics prescriptions to model the 
molecular gas phase and its CO content. The simulations 
were used to examine the CO emission properties of high-z 
mergers and starbursts (Narayanan et al. 2008c, 2009), ac¬ 
tive galactic nuclei (AGN) and quasars (Narayanan et al. 
2006, 2008b, a), as well as moderately star-forming galaxies 
(Greve & Sommer-Larsen 2008; Narayanan et al. 2011; La¬ 
gos et al. 2012; Popping et al. 2014). For example, simulating 
a large set of galaxies, Narayanan & Krumholz (2014) pro¬ 
posed a calibration scheme in which the global CO SLED 
is parameterised by the star formation rate (SFR) surface 
density, the latter being a more readily observed quantity. 
Simulations of CO have also been used to examined how well 
CO can probe the relationship between gas surface density 
and SFR, i.e., the so-called Kennicutt-Scmidt relations, and 
the validity of the CO-to-FD conversion factor (Xco) used 
to infer H 2 column densities from measured CO J = 1 — 0 
line intensities (e.g., Narayanan et al. 2011; Narayanan & 
Hopkins 2013). 

Here we present a new numerical framework for simulat¬ 
ing the molecular line emission from star-forming galaxies. 
The code, Simulator of GAlaxy Millimetre/submillimetre 
Emission (sigame), combines SPH (or grid-based) simula¬ 
tions of galaxies with sub-grid physics prescriptions for the 
H 2 /H 1 fraction and the ISM thermal balance on scales of in¬ 


dividual Giant Molecular Clouds (GMCs). sIgame accounts 
for a FUV and cosmic ray (CR) intensity field that vary with 
local SFR density within the galaxy. In this paper, we apply 
SIGAME to cosmological SPH simulations of three massive, 
normal star-forming galaxies at z = 2 (i.e., so-called main- 
sequence galaxies), and model their CO rotational line spec¬ 
trum using a publicly available 3D radiative transfer code. 
We show that SIGAME reproduces observed low-J CO line 
luminosities and provides new estimates of the aco factor 
for main-sequence galaxies at z ~ 2, while at the same time 
predicting their CO line luminosities at high-J (J up > 6) 
transitions where observations are yet to be made. 

The structure of the paper is as follows. Section 2 de¬ 
scribes the cosmological SPH simulations used, along with 
the basic properties of the three z = 2 star-forming galax¬ 
ies extracted from the simulations. A detailed description of 
SIGAME and its application to the three simulations is pre¬ 
sented in Section 3, along with a bench test of the code on 
MW-like galaxies at z = 0. The resulting CO emission maps 
and SLEDs of the z = 2 simulations are presented in Sec¬ 
tion 4, and compared to CO observations of similar galaxies 
at z ~ 1 — 2.5. In section 5, we explore how the results 
in Section 4 change when adopting alternative ISM models. 
Section 6 discusses the strengths and weaknesses of SIGAME 
in the context of other molecular line (CO) simulations and 
recent observations. Finally, in Section 7 we summarise the 
main steps of SIGAME, and list our main findings regarding 
the CO line emission from massive, star-forming galaxies at 
z ~ 2. We adopt a flat cold dark matter (ACDM) scenario 
with = 0.3, 11 a = 0.7 and h = 0.65 (Planck Collabora¬ 
tion XVI 2014). 

2 COSMOLOGICAL SIMULATIONS 
2.1 SPH simulations 

We employ a cosmological TreeSPH code for simulating 
galaxy formation and evolution, though in principle, grid- 
based hydrodynamic simulations could be incorporated 
equally well. The TreeSPH code used for the simulations is in 
most respects similar to the one described in Sommer-Larsen 
et al. (2005) and Romeo et al. (2006). A cosmic baryon frac¬ 
tion of /b = 0.15 is assumed, and simulations are initiated at 
a redshift of z = 39. The implementation of star formation 
and stellar feedback, however, has been manifestly changed. 

Star formation is assumed to take place in cold gas 
(Tk < 10 4 K) at densities rin > 1cm -3 . The star formation 
efficiency (or probability that a gas particle will form stars) 
is formally set to 0.1, but is, due to effects of self-regulation, 
considerably lower. Star formation takes place in a stochas¬ 
tic way, and in a star formation event, 1 SPH gas particle is 
converted completely into 1 stellar SPH particle, represent¬ 
ing the instantaneous birth of a population of stars accord¬ 
ing to a Chabrier (2003) stellar initial mass function (IMF; 
Chabrier 2003) - see further below. 

The implementation of stellar feedback is based on a 
sub-grid super-wind model, somewhat similar to the ‘high- 
feedback’ models by Stinson et al. (2006). These models, 
though, build on a supernova blast-wave approach rather 
than super-wind models. Both types of models invoke a 
Chabrier (2003) IMF, which is somewhat more top-heavy 
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in terms of energy and heavy-element feedback than, e.g., 
the standard Salpeter IMF. The present models result in 
galaxies characterised by reasonable z = 0 cold gas frac¬ 
tions, abundances and circum-galactic medium abundance 
properties. They also improve considerably on the “angu¬ 
lar momentum problem” relative to the models presented 
in, e.g., Sommer-Larsen et al. (2003). The models will be 
described in detail in a forthcoming paper. 


2.2 The model galaxies 

Three model galaxies, hereafter referred to as Gl, G2 and G3 
in order of increasing SFR, were extracted from the above 
SPH simulation and re-simulated using the ‘zoom-in’ tech¬ 
nique described in (e.g., Sommer-Larsen et al. 2003). The 
emphasis in this paper is on massive (M* > 5 x 10 lo M@) 
galaxies, and the three galaxies analysed are therefore larger, 
rescaled versions of galaxies formed in the 10 //iMpc cosmo¬ 
logical simulation described in Sommer-Larsen et al. (2003). 
The linear scale-factor is of the order 1.5, and since the 
CDM power spectrum is fairly constant over this limited 
mass range the rescaling is a reasonable approximation. 

Galaxy Gl was simulated at fairly high resolution, us¬ 
ing a total of 1.2 x 10 6 SPH and dark matter particles, 
while about 9 x 10 5 and 1.1 x 10 6 particles were used in 
the simulations of G2 and G3, respectively. For the Gl sim¬ 
ulation, the masses of individual SPH gas, stellar and dark 
matter particles are mspH = m* « 6.3 x 10’ and 

moM=3.5 x 10 6 A^Mq, respectively. Gravitational (cubic 
spline) softening lengths of 310, 310 and 560 h~ 1 pc, re¬ 
spectively, were employed. Minimum gas smoothing lengths 
were about 50h _1 pc. For the lower resolution simulations 
of galaxies G2 and G3, the corresponding particle masses 
are mspH = m* rts 4.7 x 10 6 /i _ 1 Mq and uidm = 2.6 x 10' 
ft _ 1 M 0 , respectively, and the gravitational softening lengths 
were 610, 610 and 1090 h - 1 pc. Minimum gas smoothing 
lengths were about 100 h - 1 pc. 

Due to effects of gravitational softening, typical ve¬ 
locities in the innermost parts of the galaxies (typically 
at radii less than about 2csph, where csph is the SPH 
and star particle gravitational softening length) are some¬ 
what below dynamical values (see, e.g. Sommer-Larsen et al. 
1998). The dynamical velocities will be of the order Vd yn = 
y/GM{R)/R, where G is the gravitational constant, R is 
the radial distance from the centre of the galaxy and M(R ) 
is the total mass located inside of R. Indeed, it turns out 
that for the simulated galaxies considered in this paper SPH 
particle velocities inside of 2esPH are only about 60-70% of 
what should be expected from dynamics. To coarsely correct 
for this adverse numerical effect, for SPH particles inside of 
2esPH the velocities are corrected as follows: For SPH par¬ 
ticles of total velocity less than Vd yn , the tangential compo¬ 
nent of the velocity is increased such that the total velocity 
becomes equal to Vd yn . Only the tangential component is in¬ 
creased in order not to create spurious signatures of merging. 
With this correction implemented, the average ratio of to¬ 
tal space velocity to dynamical velocity of all SPH particles 
inside of 2csph equals unity. 


Figure 1 shows surface density maps of the SPH gas in 
Gl, G2, and G3, i.e., prior to any post-processing by SIGAME. 
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Figure 1. SPH gas surface density maps of the three model galax¬ 
ies Gl (top), G2 (middle), and G3 (bottom) viewed face-on. The 
stellar masses and SFRs of each galaxy are indicated (see also 
Table 1). The maps have been rendered with the visualization 
tool SPLASH version 2.4.0 (Price 2007) using the gas smoothing 
lengths provided by the simulations. 
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Table 1. Physical properties of the three simulated galaxies Gl, 
G2, and G3 



SFR 

[M© yr" 1 ] 

M* 

[10 11 M 0 ] 

M S ph 

[1O 1o M 0 ] 

/sph 

Z' 

-Rcut 

[kpc] 

Gl 

40 

0.53 

2.07 

28% 

1.16 

20 

G2 

80 

1.49 

2.63 

15% 

1.97 

15 

G3 

142 

2.11 

4.66 

18% 

1.36 

20 


Notes. All quantities are determined within a radius I7 cu t, 
which is the radius where the cumulative radial stellar mass 
function of each galaxy becomes flat. The gas mass (Msph) is 
the total SPH gas mass within i? C ut- The metallicity 
( Z' = Z/Zq) is the mean of all SPH gas particles within /Gut ■ 


The gas is seen to be strongly concentrated towards the cen¬ 
tre of each galaxy and structured in spiral arms containing 
clumps of denser gas. The spiral arms reach out to a radius 
of about 20 kpc in Gl and G3, with G2 showing a more com¬ 
pact structure that does not exceed R ~ 15 kpc. Table 1 lists 
key properties of the simulated galaxies, namely their SFR, 
stellar mass (M«), SPH gas mass (Msph), SPH gas mass 
fraction (/sph = Msph/(M* + Msph)), and metallicity Z'. 
These quantities were measured within a radius (R cu t, also 
given in Table 1) corresponding to where the radial cumula¬ 
tive stellar mass function has flattened out. The metallicity 
is in units of solar metallicity and is calculated from the 
abundances of C, N, O, Mg, Si, S, Ca and Fe in the SPH 
simulations, and adjusted for the fact that not all heavy 
metals have been included according to the solar element 
abundancies measured by Asplund et al. (2009). 


The location of our three model galaxies in the SFR- 
M* diagram is shown in Figure 2 along with a sample of 
3754 1.4 < z < 2.5 main-sequence galaxies selected in near- 
IR from the NEWFIRM Medium-Band Survey (Whitaker 
et al. 2011). The latter used a Kroupa IMF but given its 
similarity with a Chabrier IMF no conversion in the stellar 
mass and SFR was made (cf., Papovich et al. (2011) and Za- 
hid et al. (2012) who use conversion factors of 1.06 and 1.13, 
respectively). Also shown is the determination of the main- 
sequence relation at 2 ~ 2 by Speagle et al. (2014), and the 
1 and 3cr scatter around it. Gl, G2, and G3 are seen to lie 
within the 3cr scatter around the z ~ 2 main sequence, albeit 
offset somewhat towards lower SFRs. This latter tendency is 
also found among a subset of CO-detected BX/BM galaxies 
at 2 ~ 2 — 2.5 (Tacconi et al. 2013), highlighted in Figure 2 
along with a handful of 2 ~ 1.5 BzK galaxies also detected 
in CO (Daddi et al. 2010). The BX/BM galaxies are selected 
by a UGR colour criteria (Adelberger et al. 2004), while the 
BzK galaxies are selected by the BzK colour criteria (Daddi 
et al. 2004). Based on the above we conclude that, in terms 
of stellar mass and SFR, our three model galaxies are repre¬ 
sentative of the star-forming galaxy population detected in 
CO at 2 ~ 2. 


3 MODELING THE ISM WITH si'game 
3.1 Methodology overview 

Here we give an overview of the major steps that go into 
SIGAME, along with a brief description of each. The full de¬ 
tails of each step are given in subsequent sections and in ap¬ 
pendices A through C. We stress, that SIGAME operates en¬ 
tirely in the post-processing stage of an SPH simulation, and 
can in principle easily be adapted to any given SPH galaxy 
simulation as long as certain basic quantities are known for 
each SPH particle in the simulation, namely: position (r), 
velocity (u), atomic hydrogen density (bh), metallicity (Z'\ 
kinetic temperature (Tk), smoothing length (h), and star for¬ 
mation rate (SFR). The key steps involved in SIGAME are: 

(i) Cooling of the SPH gas. The initially hot (Tk ~ 
10 3- ' K) SPH gas particles are cooled to temperatures typ¬ 
ical of the warm neutral medium (< 10 4 K) by atomic and 
ionic cooling lines primarily. 

(ii) Inference of the molecular gas mass fraction (fm oi = 
mH 2 /nrsPH) of each SPH particle after the initial cooling in 
step 1. /(/ 0 j for a given SPH particle is calculated by taking 
into account its temperature, metallicity, and the local CR 
and FUV radiation field impinging on it. 

(iii) Distribution of the molecular gas into GMCs. Cloud 
masses and sizes are obtained from random sampling of the 
observed CMC mass-spectrum in nearby quiescent galaxies 
and applying the local GMC mass-size relation. 

(iv) GMC thermal structure. A radial density profile is 
adopted for each GMC and used to calculate the temper¬ 
ature structure throughout individual clouds, taking into 
account heating and cooling mechanisms relevant for neu¬ 
tral and molecular gas, when exposed to the local CR and 
attenuated FUV fields. 

(v) Radiative transport of CO lines. Finally, the CO line 
spectra are calculated separately for each GMC with a ra¬ 
diative transfer code and accumulated on a common velocity 
axis for the entire galaxy. 

Determining the temperature (i) and the molecular gas 
mass fraction (ii) of a warm neutral gas SPH particle cannot 
be done independently of each other, but must be solved for 
simultaneously in an iterative fashion (see Sections 3.2 and 
3.3). As already mentioned, we shall apply SIGAME to the 
SPH simulations of galaxies Gl, G2, and G3 described in 
Section 2, and in doing so we will use them to illustrate the 
workings of the code. 


3.2 The Warm and Cold Neutral Medium 

In SPH simulations of galaxies the gas is typically not cooled 
to temperatures below a few thousand Kelvin (Springel & 
Hernquist 2003). This is illustrated in Figure 3, which shows 
the SPH gas temperature distribution (dashed histogram) in 
Gl (for clarity we have not shown the corresponding temper¬ 
ature distributions for G2 and G3, both of which are similar 
to that of Gl). Minimum SPH gas temperatures in Gl, G2 
and G3 are about 1200 K, 3100 K and 3200 K, respectively, 
and while temperatures span the range ~ 10 3 ' K, the bulk 
(~ 80 — 90%) of the gas mass in all three galaxies is at 
T k < 10 5 K. 

At these temperatures the gas will be in atomic or 
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Figure 2. Position of the three model galaxies studied here (Gl, G2 and G3 with filled a red star, a green circle and a blue square 
respectively), on a SFR-M* diagram. The grey filled contours show the z ~ 2 number density of 3754 1.4 < z < 2.5 galaxies from the 
NEWFIRM Medium-Band Survey (Whitaker et al. 2011). The SFR — M* relation at z ~ 2 as determined by Speagle et al. (2014) is 
indicated by the dashed line, with the la and 3rr scatter of the relation shown by the dot-dashed and dotted lines, respectively. Also 
shown are six z ~ 1.4 — 1.6 BzK galaxies (black circles; Daddi et al. 2010) and 14 z ~ 2 — 2.5 Bx/BM galaxies (black crosses; Tacconi 
et al. 2013). The BzK galaxies are from top to bottom: BzK—12591, BzK—21000, BzK—16000, BzK—17999, BzK—4171 and BzK—2553 
(following the naming convention of Daddi et al. 2010). 


ionised form, and H atoms that attach to dust grain sur¬ 
faces via chemical bonds (chemisorbed) will evaporate from 
the grains before H 2 can be formed. H 2 can effectively only 
exist at temperatures below ~ 10 3 K, assuming a realistic 
desorption energy of 3 x 10 4 K for chemisorbed H atoms 
(see Cazaux & Spaans 2004). The first step of SIGAME is 
therefore to cool some portion of the hot SPH gas down 
to Tk ^ 10 3 K, i.e., a temperature range characteristic of a 


warm and cold neutral medium for which we can meaning¬ 
fully employ a prescription for the formation of H 2 . 

SIGAME employs the standard cooling and heating 
mechanisms pertaining to a hot, partially ionised gas. Cool¬ 
ing occurs primarily via emission lines from H, He, C, O, 
N, Ne, Mg, Si, S, Ca, and Fe in their atomic and ionised 
states, with the relative importance of these radiative cool¬ 
ing lines depending on the temperature (Wiersma et al. 
2009). In addition to these emission lines, electron recom- 
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bination with ions can cool the gas, as recombining elec¬ 
trons take away kinetic energy from the plasma, a process 
which is important at temperatures > 10 3 K (Wolfire et al. 
2003). At similar high temperatures another important cool¬ 
ing mechanism is the scattering of free electrons off other free 
ions, whereby free-free emission removes energy from the 
gas (Draine 2011). Working against the cooling is heating 
caused by cosmic rays via the expulsion of bound electrons 
from atoms or the direct kinetic transfer to free electrons 
via Coulomb interactions (Simnett & McDonald 1969). By 
ignoring heating of the gas by photo-ionization, we shall 
adopt an assumption typically used for galactic ISM (Gnat 
& Sternberg 2007; Smith et al. 2008), though see Wiersma 
et al. (2009) who demonstrate how cooling rates are reduced 
when including photoionization in the intergalactic medium 
and proto-galaxies. 

We arrive at a first estimate of the temperature of the 
neutral medium by requiring energy rate equilibrium be¬ 
tween the above mentioned heating and cooling mechanisms: 

Pcr ,Hl — -^atoms+ions + A rec + Af_f, (1) 

where A at0 ms+ions is the cooling rate due to atomic and ionic 
emission lines, A rec and Af_f are the cooling rates from re¬ 
combination processes and free-free emission as described 
above, and For, hi is the CR heating rate in atomic, partly 
ionised, gas. The detailed analytical expressions employed 
by SIGAME for these heating and cooling rates are given in 
appendix A. 

The abundances of the atoms and ions included in 
A a toms+ions either have to be calculated in a self-consistent 
manner as part of the SPH simulation or set by hand. For 
our set of galaxies, the SPH simulations follow the abun¬ 
dances of H, C, N, O, Mg, Si, S, Ca and Fe, while for the 
abundances of He and Ne, we adopt solar mass fractions 
of 0.2806 and 10~ 4 , respectively, as used in Wiersma et al. 
(2009). 

P cr.hi depends on the primary CR ionization rate 
(Ccr), a quantity that is set by the number density of super¬ 
novae since they are thought to be the main source of CRs 
(Ackermann et al. 2013). In SIGAME, this is accounted for by 
parameterizing £cr as a function of the local star formation 
rate density (SFRD) as it varies across the simulated galaxy. 
The details of this parameterization are deferred to Section 
3.4.3. 

In eq. 1, all terms but A at0 ms+ions depend on the ion¬ 
ization fraction ( x e = n e /nm ) of the gas, which in turn 
depends on the total CR ionization rate (i.e., Ccr corrected 
for secondary ionizations of H and He), the gas tempera¬ 
ture (7k), and the Hi density (rhi) (see appendix A). The 
ionization fraction is calculated taking into account the ion¬ 
ization of H and He (with a procedure kindly provided by 
I. Pelupessy; see also Pelupessy (2005)). Since, 71 hi is set by 
the molecular gas mass fraction j/^), which in turn also 
depends on 71 (see Section 3.3 on how /(Coi is calculated), 
eq. 1 has to be solved in an iterative fashion until consistent 
values for 71, n e , and /^ ol are reached. Example solutions 
are given in Figure Al in Appendix A. 

The temperature distribution that results from solving 
eq. 1 for every SPH particle in the G1 simulation is shown in 
Figure 3 (very similar distributions are obtained for G2 and 
G3). The gas has been cooled to 71 < 10 4 K, with tempera- 



log(T k [K]) 


Figure 3. The distributions of gas kinetic temperature before 
(dashed histogram) and after (solid histogram) applying the heat¬ 
ing and cooling mechanisms of eq. 1 to galaxy Gl. The original 
hot SPH gas is seen to span a temperature range from about 
10 3 K up to ~ 10 7 K, while once the gas has been cooled the 
temperature distribution only barely exceeds ~ 10 4 K. 

tures extending down to ~ 25 K. This new gas phase repre¬ 
sents both the warm neutral medium (WNM) and the cold 
neutral medium (CNM), and from it we derive the molecular 
gas phase. 

3.3 Hi to H 2 conversion 

For the determination of the molecular gas mass fraction as¬ 
sociated with each SPH gas particle we use the prescription 
of Pelupessy et al. (2006), inferred by equating the forma¬ 
tion rate of H 2 on dust grains with the photodissociation 
rate of H 2 by Lyman-Werner band photons, and taking into 
account the self-shielding capacity of H 2 and dust extinction. 
We ignore H 2 production in the gas phase (cf., Christensen 
et al. 2012) since only in diffuse low metallicity (< 0.1 Zq) 
gas is this thought to be the dominant formation route (Nor¬ 
man & Spaans 1997), and so should not be relevant in our 
model galaxies that have mean metallicities Z' > 1 (see 
Tablet) and very little gas with Z' < 0.1 (see Figure Cl 
in Appendix C). We adopt a steady-state for the Hi—>H 2 
transition, meaning that we ignore any time dependence ow¬ 
ing to temporal changes in the UV field strength and/or 
disruptions of GMCs, both of which can occur on similar 
time scales as the H 2 formation. This has been shown to 
be a reasonable assumption for environments with metallic¬ 
ities > 0.01 Zq (Narayanan et al. 2011; Krumholz & Gnedin 
2011 ). 

The first step is to derive the FUV Held strength, Gg, 
which sets the Hi—>H2 equilibrium. In SIGAME, G 0 consists 
of a spatially varying component that scales with the local 
SFRD (SFRDi oca i) in different parts of the galaxy on top of 
a constant component set by the total stellar mass of the 
galaxy. This is motivated by Seon et al. (2011) who mea¬ 
sured the average FUV field strength in the MW (Gq iMW ) 
and found that about half comes from star light directly 
with the remainder coming from diffuse background light. 
We shall assume that in the MW the direct stellar con- 
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tribution to G' 0 MW is determined by the average SFRD 
(SFRDmw), while the diffuse component is fixed by the stel¬ 
lar mass (M*,mw). From this assumption, i.e., by calibrating 
to MW values, we derive the desired scaling relation for G' 0 
in our simulations: 


G'o — G 0 , M w 


/ SFRD local M, \ 

V ' SFRDmw ' M,, mw ) ’ 


( 2 ) 


where G^mw = 0.6 Habing (Seon et al. 2011), and M*,mw = 
6 x 10 10 Mq (McMillan 2011). For SFRDmw we adopt 
0.0024 M® yr -1 kpc -3 , inferred from the average SFR 
within the central 10 kpc of the MW (0.3 Mg yr -1 ; Hei- 
derman et al. 2010) and within a column of height equal to 
the scale height of the young stellar disk (0.2 kpc; Bovy et al. 
2012) of the MW disk. SFRDi oca i is the SFRD ascribed to a 
given SPH particle, and is calculated as the volume-averaged 
SFR of all SPH particles within a 5 kpc radius. Note, that 
the stellar mass sets a lower limit on G' 0 , which for Gl, G2, 
and G3 are 0.22, 0.62, and 0.88 Habing, respectively. 

Next, the gas upon which the FUV field impinges is as¬ 
sumed to reside in logotropic clouds, i.e., clouds with radial 
density profiles given by n(r) = riH.ext (r/R) -1 , where riH,ext 
is the density at the cloud radius R. For a logotropic density 
profile, the external density is given by riH,ext = 2/3(nu), 
where we approximate (uh) with the original SPH gas den¬ 
sity. From Pelupessy et al. (2006) we then have that the 
molecular gas mass fraction, including heavier elements than 
hydrogen, of each cloud is given by 1 : 


/mol 


m mo i Av 

- = exp -4—- 

mspH Av 


(3) 


Here (A v ) is the area-averaged visual extinction of the cloud 
and Hv tr ' is the extinction through the outer layer of neutral 
hydrogen. The area-averaged extinction, (A v ), is calculated 
from the metallicity and average cloud density, (jih): 

(A v ) = 7.21 x 10 - 22 Z'{nn)R, (4) 


where (nn)R is set by the well known density-size scaling 
relation for virialised clouds, normalised by the external 
boundary pressure, P ext , i.e.: 


(n H ) = rids 


( Pext/fcB \ 1/Z 

VlO 4 cm -3 K/ 



(5) 


For the normalization constant we adopt rids = 10 3 cm -3 , as 
inferred from studies of molecular clouds in the MW (Larson 
1981; Wolfire et al. 2003; Heyer & Brunt 2004), although 
we note that Pelupessy et al. (2006) uses 1520cm -2 . The 
external hydrostatic pressure for a rotating disk of gas and 
stars is calculated at mid-plane following Swinbank et al. 
( 2011 ): 


P tot 




fe) E 


( 6 ) 


where cr gas j_ and <t*_l are the local vertical velocity disper¬ 
sions of gas and stars respectively, and E denotes surface 


1 We will use lower case m when dealing with individual SPH 
particles. Furthermore, / alol is not to be confused with / mo i de¬ 
scribing the total molecular gas mass fraction of a galaxy and to 
be introduced in Section 4. 



Figure 4. The distribution of the H 2 gas mass fraction of the SPH 
particles in Gl calculated according to eq. 3 (solid line histogram). 
Similar distributions are found for G2 and G3. The lower limit 
on /( n(j |, defined as described in Section 3.4.1, is indicated by the 
dashed vertical line. 


densities of the same. These quantities are all calculated di¬ 
rectly from the simulation output, using the neighbouring 
SPH particles within 5 kpc, weighted by mass, density and 
the cubic spline kernel (see also Monaghan 2005). The ex¬ 
ternal cloud pressure that enters in eq. 5, is assumed to be 
equal to Ptot/(l + Qo + /3o) for relative cosmic and mag¬ 
netic pressure contributions of ao = 0.4 and /3o = 0.25 
(Elmegreen 1989; Swinbank et al. 2011). For the MW, 
Pext/fcB ~ 10 4 cm -3 K (Elmegreen 1989), but as shown in 
Figured in Appendix C, our model galaxies span a wide 
range in P ex t/fcB of ~ 10 2 — 10' cm -3 K. 

For Ai tr \ the following expression is provided by Pelu¬ 
pessy et al. (2006): 


1 (tv) 


1.086id; F u V x In 


go / Sfuv 

^S H (T k ) V Z'T k 

UH,ext/135 cm -d 


(7) 


where £fuv is the ratio between dust FUV absorption cross 
section (a) and the effective grain surface area (ad), and 
is set to (fuv = a/ao = 3. Sn = (1 + O.OlTk) -2 is the 
probability that Hi atoms stick to grain surfaces, where 
Tfc is the kinetic gas temperature (determined in an iter¬ 
ative way as explained in Section 3.2). Furthermore, v = 
n.H,ext Ra( \ T n.H,ext Ra ) 1 , and /x (set to 3.5 as suggested by 
Pelupessy et al. 2006) is a parameter which incorporates the 
uncertainties associated with primarily Sn and <Td- 

Following the above prescription SIGAME determines the 
molecular gas mass fraction and thus the molecular gas mass 
(m mo i = /moi m SPH) associated with each SPH particle. De¬ 
pending on the environment (e.g., Go, Z ', 21, P ex t,---), the 
fraction of the SPH gas converted into H 2 can take on any 
value between 0 and 1, as seen from the f' moi distribution of 
Gl in Figure 4. Overall, the total mass fraction of the SPH 
gas in Gl, G2, and G3 (i.e., within R cu t) that is converted 
to molecular gas is 29%, 52% and 34%, respectively. 
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3.4 Structure of the molecular gas 

Having determined the molecular gas mass fractions, SIGAME 
proceeds by distributing the molecular gas into GMCs, and 
calculates their masses and sizes, along with internal density 
and temperature structures, as described in the following. 


3-4-1 GMC masses and sizes 


The molecular gas associated with a given SPH particle is 
divided into GMCs by randomly sampling a power-law mass 
spectrum of the form: 

dN . . 

- -ocm G MC • (8) 

amGMC 

For GMCs in the MW disk and Local Group galaxies 0 ~ 1.8 
(Blitz et al. 2007), and is the value adopted by SIGAME in this 
paper unless otherwise stated. Lower and upper mass cut¬ 
offs at 10 4 Mg and 10 6 Mg, respectively, are enforced in or¬ 
der to span the mass range observed by Blitz et al. (2007). A 
similar approach was adopted by Narayanan et al. (2008b, c). 
Note that, in the case of G1 the upper cut-off on the GMC 
masses is in fact set by the mass resolution of the SPH simu¬ 
lation (6.3 x 10 s h _1 Mg). For Gl, typically < 30 GMCs are 
created in this way per SPH particle, while for G2 and G3, 
which were run with SPH gas particle masses almost an or¬ 
der of magnitude higher, as much as ~ 100 GMCs could be 
extracted from a given SPH particle. Figure 5 shows the re¬ 
sulting mass distribution of all the GMCs in Gl, along with 
the distribution of molecular gas mass associated with the 
SPH gas particles prior to it being divided into GMCs. The 
net effect of re-distributing the H 2 mass into GMCs is a mass 
distribution dominated by relatively low cloud masses, which 
is in contrast to the relatively flat SPH H 2 mass distribution. 
Note, the lower cut-off at ?71gmc = 10 4 Mg implies that if 
the molecular gas mass associated with an SPH particle (i.e., 
Wmoi = /moi^-SPH) is less than this lower limit it will not 
be re-distributed into GMCs. Since mspH is constant in our 
simulations (6.3 x 10 5 h -1 Mg for Gl and 4.7 x 10® Mg 
for G2 and G3) the lower limit imposed on mGMC translates 
directly into a lower limit on f mol (0.016 for Gl and 0.002 for 
G2 and G3, shown as a dashed vertical line for Gl in Figure 
4). As a consequence, 0.2, 0.005 and 0.01% of the molecu¬ 
lar gas in Gl, G2, and G3, respectively, does not end up in 
GMCs. These are negligible fractions and the molecular gas 
they represent can therefore be safely ignored. 

The GMC sizes are derived from the virial theorem, 
which relates the radius of a GMC (Rgmc) to its mass and 
external pressure (E ex t) according to: 

-Rgmc _ ( Eext/&B \ 1|/4 ( uigmc \ 

pc Vl0 4 cm- 3 K7 \290 M 0 ) ' U 


Similarly, the internal velocity dispersion ( a v ) of the clouds 
is given by: 


_ | , / 71x1 / 1 '\\ 

kms -1 ' \ 10 4 cm -3 I< 


1/4 


( Rgmc 
V P c 


1/2 


( 10 ) 


(e.g., Elmegreen 1989; Swinbank et al. 2011). Figure 6 shows 
the resulting distribution of GMC radii in Gl (solid his¬ 
togram). The minimum and maximal cloud radii found in 
Gl are ~ 0.07 pc and ~ 100 pc, respectively, and are set by 
the pressure and the imposed limits on ttigmc- 



Figure 5. The distribution of GMC masses in Gl obtained by 
applying eq. 8, with 0 = 1.8 (solid histogram) and 0 = 1.5 (dash- 
dotted histogram), compared to the total molecular gas masses 
associated with SPH particles (distribution shown as dotted his¬ 
togram). The dashed vertical line indicates the SPH gas mass 
resolution of the simulation (= 6.3 X 10 5 h 1 M 0 in the case of 
Gl). 


Observations have indicated that the shape of the GMC 
mass spectrum might be different in gas-rich systems where 
a high-pressure ISM leads to the characteristic mass of star¬ 
forming clumps of molecular gas being much higher than 
what is observed in normal spirals (e.g. Swinbank et al. 2011; 
Leroy et al. 2015). In Section 5 we therefore examine the ef¬ 
fects of adopting a more top-heavy GMC mass distribution, 
corresponding to 0 = 1.5 (shown as dot-dashed histogram in 
Figure 5). The total amount of molecular gas in our galaxies 
does not change significantly between 0 = 1.8 and 0 = 1.5, 
and the CO simulation results are robust against (reason¬ 
able) changes in the mass-spectrum. 

The GMCs are placed randomly around the position of 
their ‘parent’ SPH particle, albeit with an inverse propor¬ 
tionality between radial displacement and mass of the GMC. 
The latter is done in order to retain the mass distribution 
of the original galaxy simulation as best as possible. The 
GMCs are assigned the same bulk velocity, v, Z ', G' 0 and 
£cr as their ‘parent’ SPH particle. 


3-4-2 GMC density structure 


In order to ensure a finite central density in our GMCs, 
SIGAME adopts a Plummer radial density profile (e.g., Plum¬ 
mer 1911; Whitworth & Ward-Thompson 2001): 


ft-H 2 (R) 
cm -3 


3.55 


f m-GMc) 

V^T ) 





(ii) 


where R v is the so-called Plummer radius, which we set to 
R p = O.IRgmc- The latter allows for a broad range in gas 
densities throughout the clouds, from ~ 10 5 cm -3 in the 
central parts to a few 10 cm -3 further out. Note, eq. 11 
accounts for the presence of helium in the GMC. 

In Section 5 we examine the impact on our simulation 
results if a steeper GMC density profile with an exponent of 
—7/2 is adopted instead of the default —5/2 in eq. 11. This 
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Figure 6. The distribution of GMC radii in G1 for the adopted 
cloud mass spectrum with ft = 1.8 (solid histogram) and for a 
slightly modified spectrum with fi = 1.5 (dash-dotted histogram). 
Similar distributions are found for G2 and G3. 



^/^GMC 


Figure 7. Plummer density profiles (eq. 11) with exponents —5/2 
(solid line) and —7/2 (dashed line) for a GMC mass of 10 4 Mg 
and an external pressure of P ex t/&B = 10 4 Kcm -3 . The corre¬ 
sponding cloud radius is Pgmc = 5.9 pc (eq. 9). The vertical 
dash-dotted line marks the Plummer radius (= 0.1Pgmc)> and 
the horizontal dotted line the nn 2 value (50 cm -3 ) below which 
CO is assumed to photo-dissociate in our simulations. 


modified Plummer profile, as well as the default profile, are 
shown in Figure 7 for a GMC with mass of 10 4 M 0 and 
external pressure 10 4 Kcm -3 . 

3-4-3 GMC thermal structure 

Having established the masses, sizes, and density structure 
of the GMCs, sIgame solves for the kinetic temperature 
throughout the clouds by balancing the relevant heating and 
cooling mechanisms as a function of cloud radius. 

GMCs are predominantly heated by FUV photons (via 
the photo-electric effect) and cosmic rays. For the strength 
of the FUV field impinging on the GMCs, we take the previ¬ 


ously calculated values at the SPH gas particle position (eq. 
2). We shall assume that the CR ionization rate scales in a 
similar manner, since CRs have their origin in supernovae 
and are therefore related to the star formation and the total 
stellar mass of the galaxy: 


Ccr = Ccr.mw 


/ SFRD local M, \ 

V ' SFRDmw ’ M,, mw ) 


( 12 ) 


where SFRDi oca i, SFRDmw and M*,mw are as described in 
Section 3.3, and Ccr is scaled to the ‘canonical’ MW value of 
Ccr ,mw = 3 x 10“ 17 s _1 (e.g. Webber 1998). While the FUV 
radiation is attenuated by dust and therefore does not heat 
the GMC centres significantly, cosmic rays can penetrate 
these dense regions and heat the gas there (Papadopoulos 
et al. 2011). For this reason SIGAME attenuates the FUV field 
throughout the GMCs, while the CR ionization rate remains 
constant for a given cloud. The extinction of the FUV field 
at a certain point within a GMC is derived by integrating the 
Plummer density profile from that point and out to Rgmc- 
This H column density is converted into a visual extinction 
{Ay = Vh/2.2 x 10 21 cm -2 ) so that the attenuated FUV 
field becomes: 


G(), a tt 


= G' 0 e 


— 1.8Ay 


(13) 


where the factor of 1.8 is the adopted conversion from visual 
to FUV extinction (Black & Dalgarno 1977). 

SIGAME calculates the gas kinetic temperature through¬ 
out each GMC via the following heating and cooling rate 
balance: 


r pe + Fcr,h 2 = Ah 2 + Ago + Aqi + Aqu + A gaa -dust. (14) 


Ppe is the photo-electric heating by FUV photons, and 
Rcr,h 2 is the cosmic ray heating in molecular gas. Ah 2 is 
the cooling rate of the two lowest H 2 rotational lines (S(0) 
and S( 1)), and Ago is the cooling rate of the combined CO 
rotational ladder. A gaa _dust is the cooling rate due to in¬ 
teractions between gas molecules and dust particles, which 
only becomes important at densities above 10 4 cm -3 (e.g., 
Goldsmith 2001; Glover & Clark 2012). Aqu and Aoi are the 
cooling rates due to [Cii] and [Ol] line emission, respectively. 
The abundances of carbon and oxygen used in the prescrip¬ 
tions for their cooling rates, scale with the cloud metallicity, 
while the CO cooling scales with the relative CO to neutral 
carbon abundance ratio, set by the molecular density (see 
Appendix B). 

A ga s-duat depends on the temperature difference be¬ 
tween the gas and the dust (see eq. B8). The dust tempera¬ 
ture, Td us t, is set by the equilibrium between the absorption 
of FUV and the emission of IR radiation by the dust grains. 
We adopt the approximation given by Tielens (2005): 


Tduat „ (( Go,,tt \°' 2 

K y 1/xm ) y 10 4 Habing ) ’ 


(15) 


where Go jatt is the dust-attenuated FUV field (eq. 13) and a 
is the grain size, which we set to 1 /im for simplicity. Values 
for Tdust using eq. 15 range from 0 to 8.9 K, but we enforce 
a lower limit on Tdust equal to the z = 2 CMB temperature 
of 8.175 K. Tdu S t is therefore essentially constant (~ 8 —9K) 
throughout the inner region of the GMC models, similar 
to the value of Td us t = 8K adopted by Papadopoulos & 
Thi (2013) for CR-dominated cores. Analytical expressions 
for all of the above heating and cooling rates are given in 
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Appendix B, which also shows their relative strengths as a 
function of density for two example GMCs (Figure Bl). 

Figure C2 in Appendix C shows the resulting Tk ver¬ 
sus n h 2 behaviour for 80 GMCs spanning a broad range of 
GMC masses, metallicities, and star formation rate densi¬ 
ties. As seen in the displayed GMC models, some general 
trends can be inferred from the Tk — 71 h 2 diagrams. At fixed 
metallicity and mass, an increase in G' 0 (and therefore also 
in Ccr), leads to higher temperatures throughout the mod¬ 
els. In the outer regions this is due primarily to the increased 
photoelectric heating, while in the inner regions, heating by 
the unattenuated cosmic rays takes over as the dominat¬ 
ing heating mechanism (Figure Bl). Keeping G' 0 (and Ccr) 
fixed, lower Tk-levels and shallower Tk — uh 2 gradients are 
found in GMCs with higher metallicities. Both these trends 
are explained by the fact that the [Cn] and [Ol] cooling rates 
scale linearly with Z' (see Appendix B). 

Moving from the outskirts and inward towards the 
GMC centres, Tk drops as the attenuation of G' 0 reduces the 
photoelectric heating. However, the transition from cooling 
via [Cn] to the less efficient cooling mechanism by CO lines, 
causes a local increase in Tk at riH 2 ~ 10 3 — 10 4 ’ 5 cm -3 , 
as also seen in the detailed GMC simulations by Glover & 
Clark (2012). The exact density at which this ‘bump’ in Tk 
occurs depends strongly on the mass of the GMC. Our choice 
of the Plummer model for the radial density profile means 
that the extinction, Ay, at a certain density increases with 
GMC mass. This in turn decreases the FUV heating, and as 
a result the Tk — «h 2 curve moves to lower densities with 
increasing GMC mass. 

As the density increases towards the cloud centres (i.e., 
tih 2 > 10 4 - 5 cm -3 ) molecular line cooling and also gas-dust 
interactions become increasingly efficient and start to dom¬ 
inate the cooling budget. The Tk — n h 2 curves are seen to 
be insensitive to changes in Z' , which is expected since the 
dominant heating and cooling mechanisms in these regions 
do not depend on Z'. Eventually, in the very central regions 
of the clouds, the gas reaches temperatures close to that of 
the ambient CMB radiation field, irrespective of the overall 
GMC properties and the conditions at the surface (Figure 
Bl). 

3.4-4 GMC grid models 

The Tk — n h 2 curve for a given GMC is determined by the 
following quantities: 

• Go and Ccr, which govern the gas heating and are set 
by the local star formation rate density (SFRDi oca i) and the 
total stellar mass (M„) according to eqs. 2 and 12. 

• tugmc and P ex t which determine the effective radius of 
a cloud (eq. 9) and thus its density profile (eq. 11). 

• The local metallicity (Z 1 ), which influences the fraction 
of H 2 gas and plays an important role in cooling the gas. 

These local parameters together with the most impor¬ 
tant global parameters used by SIGAME are listed in Ta¬ 
ble 2. The GMC ensemble distributions of G' 0 , Ccr, uigmc, 
Text/^B, and Z' for each of the galaxies Gl, G2 and G3 are 
shown in Figure Cl. 

There are more than 100,000 GMCs in a single model 
galaxy and, as Figure Cl shows, they span a wide range in 


G 0 , Ccr, w!gmc, Text/te, and Z'. Thus, in order to shorten 
the computing time, we calculated Tk — tih 2 curves for a set 
of 630 GMCs, chosen to appropriately sample the distribu¬ 
tions at certain grid values (listed in Table 3, and marked by 
vertical black lines in Figure Cl). Every GMC in our simu¬ 
lations was subsequently assigned the Tk — nn 2 curve of the 
GMC grid model closest to it in the (Go, w.gmc, Z') parame¬ 
ter space. In the default GMC grid, we keep Text/te fixed to 
a MW-like value of 10 4 Kcm -3 , thereby also anchoring the 
scaling relations in eq. 9 and 10 for size and velocity disper¬ 
sion. This is done to minimise the number of radiative trans¬ 
fer calculations, effectively keeping the amount of molecular 
gas mass dependent on local pressure, but removing local 
pressure as a free parameter in, and hence simplifying, the 
calculation of CO emission. 

In addition to the default grid described above, we made 
two separate tests to explore the GMC parameter space 
more fully. These tests are described in the following and 
results of their application are presented in Section 5. The 
external cloud pressure, P e xt/feB, spans a wide range from 
~ 10 3 to 10 9 Kcm -3 as shown in Figure Cl. In a separate 
test, we therefore constructed the same GMC grid, but for 
fixed pressures of Pext/fcB = 10 5 ' 5 and 10 6 ' 5 Kcm -3 and 
interpolated among the resulting three values of P ex t/kB 
([10 4 ,10 5 ' 5 ,10 6 ' 5 ] Kcm -3 ) as a demonstration of how to in¬ 
corporate local pressure in SIGAME. We also test the ef¬ 
fect of adopting the alternative radial GMC density profile 
defined in Section 3.4.2. Tk — tih 2 curves were calculated 
for each possible combination of the parameter grid values 
listed in Table 3 for both density profiles, giving us a total 
of 3 x 2 x 630 = 3780 GMC grid models. 


3.5 Radiative transfer of CO lines 

SIGAME assumes a fixed CO abundance equal to the Galac¬ 
tic value of [CO/H 2 ]= 2 x 10 -4 (Lee et al. 1996; Sofia et al. 
2004, and see Section 6 for a justification of this value) ev¬ 
erywhere in the GMCs except for riH 2 < 50 cm -3 , where CO 
is not expected to survive photo-dissociation processes (e.g. 
Narayanan et al. 2008c). SIGAME calculates the CO line ra¬ 
diative transfer for each GMC individually and derives the 
CO line emission from the entire galaxy. 

3.5.1 Individual GMCs 

For the CO radiative transfer calculations we use a slightly 
modified version of the Line Modeling Engine (lime ver. 
1.4; Brinch & Hogerheijde 2010) - a 3D molecular excita¬ 
tion and radiative transfer code, lime has been modified in 
order to take into account the redshift dependence of the 
CMB temperature, which is used as boundary condition for 
the radiation field during photon transport, and we have also 
introduced a redshift correction in the calculation of phys¬ 
ical sizes as a function of distance. We use collision rates 
for CO (assuming H 2 is the main collision partner) from 
Yang et al. (2010). In lime, photons are propagated along 
grid lines defined by Delaunay triangulation around a set 
of appropriately chosen sample points in the gas, each of 
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Table 2. Parameters, global and local, used by SIGAME, together with relevant equations in this paper. 


Global parameters [CO/H 2 ], ft (GMC mass spectrum), GMC density profile Eqs. 8 and 11 

Local parameters mcMC) Gq, Ccr? Z ', P ex t, SFRD Eqs. 8, 2, 2 and 6 

Derived internal GMC parameters njj, 7k, x e , cr v Eqs. 11 and 14 


Table 3. Grid parameter values 


G' 0 [Habing] 0.5, 1, 4, 7, 10, 13, 16, 19, 23, 27 

log(m G MC [M©]) 4.0, 4.25, 4.5, 4.75, 5.0, 5.25, 5.5, 5.75, 6.0 

log (Z/Zq) -1, -0.5, 0, 0.5, 1, 1.4, 1.8 

log(Pext/fc B [K cm -3 ]) 4.0, 5.5, 6.5 


which contain information on nn 2 , Tk, a v , [CO/H 2 ] and v. 
SIGAME constructs such a set of sample points throughout 
each GMC: about 5000 points distributed randomly out to 
a radius of 50 pc, i.e., beyond the effective radius of typical 
GMCs in Gl, G2, and G3 (Figure 6) and in a density regime 
below the threshold density of 50 cm -3 adopted for CO sur¬ 
vival (see previous paragraph). The concentration of sample 
points is set to increase towards the centre of each GMC 
where the density and temperature vary more drastically. 

For each GMC, lime generates a CO line data cube, 
i.e., a series of CO intensity maps as a function of velocity. 
The velocity-axis consists of 50 channels, each with a spec¬ 
tral resolution of 1.0 km s - , thus covering the velocity range 
v = [—25,+25] kins -1 . The maps are 100pc on a side and 
split into 200 pixels, corresponding to a linear resolution of 
0.5pc/pixel (or an angular resolution of 5.9 x 10 -5 "/pixel 
at z = 2). Intensities are corrected for arrival time delay and 
redshifting of photons. 

Figure C3 shows the area- and velocity-integrated CO 
Spectral Line Energy Distributions (SLEDs) for the same 
80 GMCs used in Section 3.4.3 to highlight the Tk — riH 2 
profiles (Figure C2). The first thing to note is that the CO 
line fluxes increase with jit-gmc , which is due to the increase 
in size, i.e., surface area of the emitting gas, with cloud mass 
(Section 3.4.1). Turning to the shape of the CO SLEDs, a 
stronger G' 0 (and £cr) increases the gas temperature and 
thus drives the SLEDs to peak at higher J-transitions. Only 
the higher, J up > 4, transitions are also affected by metallic- 
ity, displaying increased flux with increased Z'. For GMCs 
with high Go, the high metallicity levels thus cause the CO 
SLED to peak at J up > 8. 


3.5.2 The effects of dust 

Dust absorbs the UV light from young O and B stars and 
re-emits in the far-IR, leading to possible TR pumping’ of 
molecular infrared sources. However, due to the large vibra¬ 
tional level spacing of CO, the molecular gas has to be at a 
temperature of at least 159 K, for significant IR pumping of 
the CO rotational lines to take place, when assuming a max¬ 
imum filling factor of 1, as shown by Carroll & Goldsmith 
(1981). Most of the gas in our GMC models is at temper¬ 
atures below 100 K, with only a small fraction of the gas, 
in the very outskirts, of the GMCs reaching Tk > 159 K, 
as seen in Figure C2 in Appendix C. This happens only if 


the metallicity is low ( Z' < 0.1) or in case of a combination 
between high FUV field ( G' 0 > 4) and moderate metallicity. 

In principle, the high-J CO lines could be subject to 
extinction by dust, but this effect is significant only in 
extremely dust-enshrouded sources (Papadopoulos et al. 
2010), and therefore unlikely to be relevant for our simu¬ 
lations. 

While lime is capable of including dust in the radiative 
transfer calculations, provided that a table of dust opaci¬ 
ties as function of wavelength be supplied together with the 
input model, we have chosen not to include dust in the sim¬ 
ulations presented here. 

3.5.3 CO emission maps 

A given GMC is assigned the CO emission line profile of 
the GMC model that corresponds to the nearest grid point 
in the (G' 0 , m gmc, Z') parameter space. A spatial grid of 
400 x 400 pixels is overlayed on each galaxy viewed face-on, 
and the line profiles of all GMCs within each pixel are added 
to a common velocity axis, thus resulting in a position- 
velocity datacube. CO moment 0 maps are then constructed 
by integrating the flux within each pixel in velocity. These 
maps are 40 kpc on each side and thus have a resolution of 
~ lOOpc/pixel. By expanding the pixel size to encompass 
the entire galaxy, the global CO SLED is derived. 

The approach described above assumes that the molec¬ 
ular line emission from each GMC is radiatively decoupled 
from all other GMCs, due to the significant velocity gradi¬ 
ents across the galaxies (~ 200 — 500kms -1 ) and the rel¬ 
atively small internal line widths of the individual GMCs 
(~ 2 — 34kms -1 for the GMCs in Gl, G2 and G3). 


4 SIMULATING MASSIVE Z = 2 MAIN 
SEQUENCE GALAXIES 

In this section we examine the H 2 surface density and CO 
emission maps resulting from applying SIGAME to the three 
SPH galaxy simulations Gl, G2, and G3 at z = 2 (Section 
2), and we compare with existing CO observations of main 
sequence galaxies at z ~ 1 — 2.5. We also tested sigame on 
three MW-like galaxies (see Appendix D), finding a general 
agreement with the CO line observations of MW and other 
similar local galaxies. As mentioned in previous sections, our 
default grid will be that corresponding to a Plummer density 
profile and a pressure of Pext/ks = 10 4 cm -3 K, combined 
with a GMC mass spectrum of slope [3 = 1.8, unless other¬ 
wise stated. 

4.1 Total molecular gas content and H 2 surface 
density maps 

The total molecular gas masses of Gl, G2 and G3 - ob¬ 
tained by summing up the GMC masses associated with all 
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Figure 8. Top row: Molecular surface density maps of our model galaxies seen face-on. The maps have been smoothed using a circular 
Gaussian with full width at half maximum (FWHM) of 3 pixels corresponding to 0.24 kpc. The molecular gas surface density maps 
are seen to trace the spiral arms and inner disk as well as a few dense clumps further out. Bottom row: azimuthally averaged radial 
profiles of the molecular (solid curve) and Hi +Hii (dashed curve) gas surface densities, of the mean metallicity (dotted curve) and of 
G' 0 (dot-dashed curve) - determined from 50 radial bins stretching from 0 to 15 kpc from the centre of each galaxy. The molecular gas 
is the result of applying the recipes in Section 3.3 to the SPH simulations presented in Section 2.2, while the Hi + Hu gas is the initial 
SPH gas mass minus the derived molecular gas mass (both include the contribution from helium). We estimate G' 0 by averaging over the 
FUV fields impinging on all GMCs in each radial bin. 


SPH particles within each galaxy (i.e., M mo i = ^tug mc = 
E /moirasPH) - are 7.1 x 10 9 , 1.7 x 10 10 , and 2.1 x 10 10 M 0 , 
respectively, corresponding to about 34, 59 and 45 % of 
the original total SPH gas masses of the galaxies within 
R cut = 20 kpc. 2 The global molecular gas mass fractions 
(i.e., /moi = M mo i/(M* + M mo i) are 11.8, 10.2 and 9.1% for 
Gl, G2, and G3, respectively. 

This is ~ 4 — 5 x below the typical molecular gas mass 
fraction (~ 40—60 %) inferred from CO observations of main 
sequence galaxies at z ~ 1 — 3 with similar stellar masses 
and star formation rates as our simulations (e.g., Daddi et al. 
2010; Magnelli et al. 2012; Tacconi et al. 2013). 

The reason for this discrepancy is the low SPH gas mass 
fractions of our simulated galaxies to begin with (/sph = 
Msph/ (M* + Msph) = 18 — 28%), which obviously restricts 
/moi to lower values. Such low gas fractions have been ob¬ 
served in local star forming galaxies. For example, assum¬ 
ing a MW-like aco factor, Saintonge et al. (2011) derived 
/moi ~ 6.7±4.5 % for a sample of 119 CO(l—0) detected nor¬ 
mal star forming galaxies and a stack of 103 non-detections. 


2 These global molecular-to-SPH gas mass fractions are calcu¬ 
lated as Mmol/MgpH = 52 m GMC / 52 m SPHi where the sums are 
over all SPH particles. 


Figure 8 (top panels) shows the molecular gas surface 
density (E mo i) maps of Gl, G2, and G3, where E mo i is cal¬ 
culated over pixels 80 pc x 80 pc in size. The pixel size was 
chosen in order to avoid resolving the GMCs, which typically 
have sizes g 40 pc (see Figure 6). The molecular gas is seen 
to extend out to radii of ~ 10 kpc and beyond, but generally 
the molecular gas concentrates within the inner regions of 
each galaxy. The distribution of molecular gas broadly fol¬ 
lows the central disk and spiral arms where the SPH surface 
density (Esph) is also the highest (Figure 1). The correspon¬ 
dence is far from one-to-one, however, as seen by the much 
larger extent of the SPH gas, i.e., regions where H 2 has not 
formed despite the presence of atomic and ionised gas. This 
point is further corroborated in the bottom panels of Figure 
8, which show azimuthally-averaged radial surface density 
profiles of the molecular gas and of the Hi + Hu gas. The 
latter is simply the initial SPH gas mass with the molecular 
gas mass subtracted, and has been corrected (just like the 
molecular gas phase) for the mass contribution from helium. 

In order to compare Ehi+hii and E mo i with observa¬ 
tions, we have set the radial bin width to 0.5 kpc - the 
typical radial bin size used for nearby spirals in the work 
of Leroy et al. (2008) 3 . The radially binned S mo i reaches 

3 The H 2 surface density maps in Figure 8 are averaged over areas 
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Figure 9. The top three rows show the moment 0 maps of the CO(l—0), CO(3—2) and CO(7—6) emission from Gl, G2 and G3. The 
CO maps have been smoothed using a circular Gaussian with full width at half maximum (FWHM) of 3 pixels corresponding to 0.24 kpc, 
and as with the SPH and molecular gas surface density maps, a logarithmic scale has been applied in order to better display extended 
emission. The bottom row shows the azimuthally averaged CO 3—2/1—0 and 7—6/1—0 brightness temperature line ratios (denoted r ^2 
and 7*76, respectively) as functions of projected radius for each of the three galaxies. A radial bin-size of 0.5 kpc was used. Also shown 
are the azimuthally averaged radial profiles of CO-to-H 2 conversion factor olqo(R) = £jj 2 /1cO(i—0) units of Mq pc -2 (Kkms -1 ) -1 
with a dashed line indicating the global aco factor and a dotted line for the disk-averaged a-co factor (see Section 4.2). 
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~ 800—1000 Mg pc -2 in the central regions of our simulated 
galaxies, which is comparable to observational estimates of 
£ mo i (of several 100 M@ pc -2 ) towards the centres of nearby 
spirals (Leroy et al. 2008). In all three galaxies, the Hl+Hn 
surface density dips within the central ~ 1 — 2 kpc, coincid¬ 
ing with a strong peak in £ mo i. Thus, despite the marked 
increase in the FUV radiation field towards the centre, the 
formation of H 2 driven by the increase in gas pressure is 
able to overcome photodestruction of H 2 through absoption 
of Lyman or Werner band radiation. The central H 2 surface 
densities are similar for all galaxies and is a direct conse¬ 
quence of very similar SPH gas surface densities in the cen¬ 
tre combined with molecular gas mass fractions approaching 
1. From R ~ 2 kpc and out to ~ 10 kpc, the Hi + Hii surface 
density remains roughly constant with values of ~ 40, ~ 70 
and ~ 100 M® pc -2 for Gl, G2 and G3, respectively. 

Radial profiles of the Hi and molecular gas surface 
density that are qualitatively very similar to our simula¬ 
tions have been observed in several nearby star-forming disk 
galaxies (e.g., Leroy et al. 2008; Bigiel et al. 2008). In local 
galaxies, however, the Hi surface density, including helium, 
rarely exceeds ~ 10 Mq pc -2 , while in our simulations we 
find Hi -)-Hii surface densities that are 4— 10 x higher, which 
is due to the substantial fraction of ionised gas in our simu¬ 
lated galaxies. 

4.2 CO line emission maps and resolved 
excitation conditions 

Moment 0 maps of the CO(l—0), CO(3—2) and CO(7— 6 ) 
emission from Gl, G2 and G3 are shown in Figure 9. Both 
the CO(l—0) and CO(3—2) emission are seen to trace the H 2 
gas distribution well (Figure 8 ), while the CO(7— 6 ) emission 
only trace gas in the central ~ 7 kpc of the galaxies. 

Also shown in Figure 9 (bottom row) are the az- 
imuthally averaged CO 3—2/1—0 and 7—6/1—0 brightness 
temperature line ratios (denoted r 32 and respectively) 
as a function of radius for Gl, G2 and G3. The profiles 
show that the gas is more excited in the central ~ 5 kpc, 
where typical values of r 32 and rye are ~ 0.55 — 0.65 and 
~ 0.02 — 0.08, respectively, compared to r 32 ~ 0.5 and 
r 76 < 0.01 further out in the disk. This radial behaviour 
of the line ratios does not reflect the H 2 gas surface density, 
which peaks towards the centre rather than flattens, and 
gradually trails off out to R ~ 12 kpc instead of dropping 
sharply at R ~ 4 —6 kpc (Figure 9). Rather, r 32 and r-jQ seem 
to mimick the radial behaviour of G' 0 (and thus £cr), which 
makes sense since Go and £cr are the most important fac¬ 
tors for the internal GMC temperature distribution (Figure 
C2). The central values for r 32 and increase when going 
from Gl to G3, as expected from the elevated levels of star 
formation density (and of G' Q and Ccr, accordingly) in their 
central regions. Beyond ~ 6 kpc the line ratios are constant 
(~ 0.5 and < 0.01) and the same for all three galaxies, due 
to relatively similar G' 0 (and Ccr) there. 

The decrease in r 32 towards the centre has also been ob¬ 
served in nearby galaxies. For M51, a spiral galaxy with a 

80 X 80 pc in size, and therefore give higher peak surface densities 
than the radial profiles which are averaged over ~ 0.5 kpc wide 
annuli. 


smaller SFR than our model galaxies (see Appendix D), Vla- 
hakis et al. (2013) measured a median ratio of r 32 = 0.54 for 
pixels covering the central kpc region, but typical ratios of 
0.2 — 0.4 in the arm and inter-arm regions. In comparison, 
our model galaxies display larger r 32 ratios and a slightly 
less pronounced drop of 0.1 or less when going from the 
central kpc region to the outskirts of the disk. Mao et al. 
(2010) observed 125 nearby galaxies of different types and 
found global r 32 values to be 0.61 ± 0.16 in normal galaxies, 
and > 0.89 in starbursts and (U)LIRGs. Iono et al. (2009) 
and Papadopoulos et al. (2012) found slightly lower r 32 in 
their samples of (U)LIRGs, with mean values of 0.48 ± 0.26 
and 0.67 ± 0.62 respectively. The r 32 profiles of our model 
galaxies reveal more highly excited gas in their centres than 
in the disks, with central values similar to the observed val¬ 
ues towards (U)LIRGs by Papadopoulos et al. (2012) but 
below those reported by Mao et al. (2010). Finally, Geach 
& Papadopoulos (2012) employed Large Velocity Gradient 
(LVG) radiative transfer models to examine r 32 in differ¬ 
ent environments. For quiescent clouds similar to very low- 
excitation gas clouds found in M31, they showed that r 32 
can be as low as ~ 0.13, while dense (n(H 2 ) > 10 4 cm -3 ) 
star-forming clouds typically have r 32 ~ 0.88. Most likely 
a result of their moderate star formation rates, the r 32 ra¬ 
dial profiles of our model galaxies lie in between these more 
extreme cases. 

4.3 The CO-to-H 2 conversion factor 

The CO-to-H 2 conversion factor (aco) connects CO(l—0) 
line luminosity (in surface brightness temperature units) 
with the molecular gas mass (M mo 1 ) as follows: 


U LU j ! ? \- LKJ J 

^CO(l-O) 

From the CO(l — 0) surface brightness and H 2 surface den¬ 
sity maps of our simulated galaxies, we calculate the average 
qco within radial bins from the galaxy centres (Figure 9). 
The resulting radial profiles show that aco is essentially 
constant as a function of radius, taking on values in the 
range ~ 1.3— 1.8 Mq pc -2 (Kkms -1 ) -1 across all three 
galaxies, except in the outskirts of G2 where aco shoots up 
to higher values. The high values of aco in the outskirts 
of G2 reflect the very low molecular gas surface densities 
(see Figure 8 ) and resulting lack of CO line emission, ren¬ 
dering aco a meaningless quantity in this region. In order 
to compare with the stellar properties in Table 1, we will use 
only the region within R = 15 kpc for the further analysis of 
CO line emission in G02. Small systematic changes in aco 
with radius are seen in all three simulated galaxies as aco is 
systematically below (above) the global aco value, marked 
with dashed lines, at <; 5kpc (^ 5kpc). Blanc et al. (2013) 
measured a drop in aco by a factor of two when going from 
R ~ 7 kpc to the R < 2 kpc central region of the Sc galaxy 
NGC 628, assuming a constant gas depletion timescale when 
converting SFR surface densities into gas masses. For com¬ 
parison, the radial aco profiles of our model galaxies only 
drop by about 10% from R> 7 kpc to the R < 2 kpc central 
region. 

Sandstrom et al. (2013) measured and examined the re¬ 
solved aco values in a sample of 26 nearby spiral galaxies. 
A disk-averaged aco value was calculated for each galaxy, 
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by tiling them with pixels of spacing 37.5” (corresponding 
to pixel sizes of 0.7 — 3.9 kpc 2 ) and calculating the mean 
of aco values for pixels with a sufficient signal-to-noise ra¬ 
tio. It was found that, on average, the galaxies have a cen¬ 
tral (R < 1 kpc) aco value a factor of two below the disk- 
averaged aco value. We calculate disk-averaged aco for our 
model galaxies in a similar way by tiling the galaxies seen 
face-on with 1 kpc 2 -sized pixels and taking the average over 
pixels within the cut-out radii R cu t given in Table 1. The 
resulting values are shown with dotted lines in Figure 9 and 
are no more than ~ 1.2x the central (R < 1 kpc) aco- Our 
simulated galaxies, therefore do not quite reproduce the drop 
in aco typically observed when going from the disk to the 
central regions of local, spiral galaxies. Although, we note 
that for three galaxies from the Sandstrom et al. (2013) sam¬ 
ple (NGC3938, NGC3077 and NGC4536) aco changes by 
only 10 — 16 % from the disk to the centre, which is in line 
with our simulations. 

The aco factor is expected to depend on Z' , as higher 
metallicity means higher C and O abundances as well as 
more dust that helps shield CO from photodestruction by 
FUV light, thereby leading to a possibly lowering of aco- 
A comparison of the aco radial profiles with those of G' 0 
and Z' in the bottom panel of Figure 8, suggests that the 
transition in aco is caused by a change in Go rather than 
a change in Z', since aco and G' 0 generally start to drop at 
around R ~ 6 kpc while Z' already drops drastically at 1 kpc 
from the centre. Our modeling therefore implies that aco 
is controlled by G' 0 rather than Z' in normal star-forming 
galaxies at z ~ 2, in agreement with the observations by 
Sandstrom et al. (2013) who do not find a strong correlation 
with Z'. 

From the total molecular gas masses and CO(l — 0) lu¬ 
minosities of Gl, G2, and G3, we derive global aco factors of 
aco = 1.6, 1.5 and 1.4M© pc -2 (Kkms -1 ) -1 , respectively. 
These values are lower (by a factor ~ 3) than the inner disk 
MW value (aco.Mw — 4.3 ± 0.1 M@ pc -2 (Kkms -1 ) -1 ), 
and closer to the typical mergers/starburst aco-values (~ 
0.2 x oco.mw) inferred from CO dynamical studies of lo¬ 
cal ULIRGs (e.g., Solomon et al. 1997; Downes & Solomon 
1998; Bryant & Scoville 1999) and 2 ~ 2 SMGs (e.g., Tac- 
coni et al. 2008). Our aco-values are also below those in¬ 
ferred from dynamical modeling of z ~ 1.5 BzKs (aco = 
3.6 ± 0.8Mg pc -2 (Kkms -1 ) -1 ; Daddi et al. 2010). 

The aco factors of the 2 ~ 1 — 1.3 star-forming galax¬ 
ies studied by Magnelli et al. (2012) occupy the same re¬ 
gion of the M*- SFR plan as our model galaxies, but have 
aco factors at least a factor ~ 6 higher. Their aco are 
~ 5—20 M© pc -2 (K km s -1 ) -1 when converting dust masses 
to gas masses using a metallicity-dependent gas-to-dust ra¬ 
tio, and Magnelli et al. (2012) further find that the same 
galaxies have relatively low dust temperatures of < 28 K 
compared to the galaxies further above the main sequence 
of Rodighiero et al. (2010). The fact that the main-sequence 
galaxies of Magnelli et al. (2012) are all of near-solar metal¬ 
licity as ours, cf. Table 1, suggests that other factors, such 
as dust temperature caused by strong FUV radition, can be 
more important than metallicity in regulating aco, in line 
with our study of the aco radial profiles. 


4.4 Global CO line luminosities and spectral line 
energy distributions 

The global CO SLEDs of Gl, G2 and G3 are shown 
in Figure 10 in three different incarnations: 1) total CO 
line luminosities {L' CO] ), 2) brightness temperature ra¬ 
tios (Tb.cOj.j.JTb.cOx.o = l ' C03 3 1 / L ' COl0 ) and 3) line 
intensity ratios (/cOj.j.J lco ll0 = L co 3 31 / L' COl 0 x 

( i 'co.j i j_i/z'cOi, 0 )”). We have also compiled relevant CO 
observations of normal star-forming galaxies at 2 ~ 1 — 3 
in order to facilitate a comparison with our simulated CO 
SLEDs. In addition to the 2 ~ 1 — 1.5 BzK and 2 ~ 2 — 2.5 
BX/BM samples by Daddi et al. (2010) and Tacconi et al. 
(2013), respectively, this includes 7 star-forming galaxies 
(SFGs) at 2 ~ 1.2 (Magnelli et al. 2012) and 39 SFGs at 
2 ~ 1 — 1.5 (Tacconi et al. 2013). The galaxies in the sample 
from Magnelli et al. (2012) are all detected in CO(2—1) as 
well as in Herschel/ PACS bands, and have log(M*/M©) = 
10.36 - 11.31 and SFR~ 29 - 74M© yr -1 . The 2 ~ 1 - 1.5 
SFG sample from Tacconi et al. (2013) comes from the Ex¬ 
tended Growth Strip International Survey (EGS), is covered 
by the CANDELS and 3D-HST programs (J — H bands and 
Ha respectively), and has log(M*/M©) = 10.40— 11.23 and 
SFR~ 28- 630 M 0 yr -1 . 

The 2 ~ 2—2.5 BX/BM galaxies of Tacconi et al. (2013) 
are not only closest in redshift to our simulated galaxies 
but also occupy the same region of the SFR—M* plane (see 
Figure 2). We find that G2 and G3, the two most massive 
(M* > 10 11 Mq) and star-forming (SFR > 80 M© yr -1 ) 
galaxies of our simulations, have CO(3—2) luminosities of 

6.4 x 10® and 9.2 x 10 9 Kkms -1 pc 2 , respectively, in the 
range of CO(3—2) luminosities of the BX/BM galaxies. Gl 
is just below this range and about a factor of 4 below the av¬ 
erage BX/BM luminosity (~ 1.2±0.9 x 10 10 Kkms -1 pc 2 ). 
All three simulated galaxies have CO(2—1) and CO(3—2) 
luminosities consistent with those of the 2 ~ 1 — 1.5 SFGs 
observed by Magnelli et al. (2012) and Tacconi et al. (2013), 
which span a range in CO(2—1) and CO(3—2) luminosity of 
(3.6- 12.5) x 10 9 and (1.4-41.3) x 10 9 K km s -1 pc 2 , re¬ 
spectively. Finally, we see that the 2 ~ 1.5 BzK galaxies in 
general have higher CO line luminosities across all observed 
CO transitions (up to J up = 5), although at the J = 2 — 1 
and J = 3 — 2 transitions there is overlap with G3. As men¬ 
tioned in Section 4.1, the molecular gas mass fractions of 
our galaxies is a factor ~ 4 — 5 below the mean of the ob¬ 
served galaxies at 2 ~ 1 — 2.5 with which we compare, and 
we propose this as the main reason for the comparatively 
low CO luminosities. 

Differences in the global CO excitation conditions be¬ 
tween Gl, G2 and G3 are best seen in the CO(l—0) nor¬ 
malised luminosity (and intensity) ratios, i.e., middle (and 
bottom) panel in Figure 10. The CO SLEDs of all three 
galaxies follow each other quite closely up to the J = 3 — 2 
transition where the SLEDs all peak; at 3 < J up < 7 the 
SLEDs gradually diverge. While the metallicity distributions 
in Gl, G2 and G3 are similar (see Figure Cl), the rise in G'o 
presents a likely cause to the increasing high-J flux when go¬ 
ing from Gl to G3, as higher Go leads to more flux primarily 
in the J > 4 transitions (see Figure C3). 

Our simulated galaxies are seen to have CO 2—1/1—0, 
3—2/1—0, and 5—4/1—0 brightness temperature ratios of 
r 2 i — 1, r "32 — 0.6, and rs 4 ~ 0.15, respectively. The first 
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Figure 10. Global CO SLEDs of our three model galaxies Gl, G2 and G3 shown as red (dashed), green (dash-dot) and blue (dash-dot- 
dot-dot) curves, respectively. The SLEDs are given as absolute line luminosities in units of Kkms -1 pc -2 (top panel), as brightness 
temperature ratios normalised to the CO(l—0) transition (middle panel), and as velocity-integrated intensity ratios normalised to 
CO(l—0) (bottom panel). The model CO SLEDs are compared with observations of z ~ 1.4—1.6 BzK galaxies (open circles; Dannerbauer 
et al. 2009; Daddi et al. 2010, 2015; Aravena et al. 2010, 2014), z ~ 1 — 1.5 star-forming galaxies (CO(3—2); empty triangles; Tacconi 
et al. 2013), z ~ 1 — 1.3 star-forming galaxies (CO(l—0) and CO(2—1); right-facing triangles; Magnelli et al. 2012), and z ~ 2 — 2.5 
BX/BM galaxies (CO(3—2); crosses; Tacconi et al. 2013). Four BzK galaxies (BzK—16000, BzK—4171, BzK—21000 and BzK—610) have 
been observed in CO(l — 0) and at least one additional transition to date, and are highlighted in the bottom two panels by connecting 
dashed lines and individual symbols. Also shown in the bottom panel, with dotted lines, are the line ratio predictions of Narayanan Sz 
Krumholz (2014) (D14), calculated for the Esfr of our galaxies, Gl to G3 from bottom to top (see Section 6). 
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two ratios compare extremely well with the line ratios mea¬ 
sured for BzK—4171 and BzK—21000, i.e., V 2 i — 0.9 — 1, 
and 7'32 — 0.5 — 0.6 (Dannerbauer et al. 2009; Aravena et al. 
2010; Daddi et al. 2010, 2015), and suggest that our simu¬ 
lations are able to emulate the typical gas excitation con¬ 
ditions responsible for the excitation of the low-J lines in 
normal 2 ~ 1 — 3 SFGs. In contrast, rs 4 = 0.3 — 0.4 ob¬ 
served in BzK—4171 and BzK—21000 (Daddi et al. 2015), 
is nearly 2 x higher than our model predictions. Daddi et al. 
(2015) argue that this is evidence for a component of denser 
and possibly warmer molecular gas, not probed by the low- 
J lines. In this picture, we would expect CO(4—3) to probe 
both the cold, low-excitation gas as well as the dense and 
possibly warm star-forming gas traced by the CO (5—4) line, 
and we would expect CO(6—5) to be arising purely from this 
more highly excited phase and thus departing even further 
from our models. 

However, significant scatter in the CO line ratios of 
main-sequence galaxies is to be expected, as demonstrated 
by the significantly lower line ratios observed towards 
BzK—16000: r 2 i — 0.4, r 32 — 0.3, and rs 4 ~ 0.1 (Aravena 
et al. 2010; Daddi et al. 2010, 2015) and, in fact, the aver¬ 
age 7*54 for all three BzK galaxies above (~ 0.2; Daddi et al. 
2015) is consistent with our models. It may be that G2 and 
G3, and perhaps even Cl, are more consistent with the CO 
SLEDs of the bulk 2 ~ 1 — 3 main-sequence galaxies. We 
stress, that to date, no z ~ 1 — 3 main sequence galaxies 
have been observed in the J = 4 — 3 nor the 6 — 5 transi¬ 
tions, and observations of these lines, along with low- and 
high-J lines in many more BzK and main-sequence 2 ~ 1 —3 
galaxies are needed in order to fully delineate the global CO 
SLEDs in a statistically robust way. 


5 TESTING DIFFERENT ISM MODELS 

In this section we investigate the effects on our simulation 
results when adopting i) a more top-heavy CMC mass spec¬ 
trum with a slope of /3 = 1.5, ii) a steeper CMC density 
profile, i.e., a Plummer profile exponent of —7/2, and iii) a 
CMC model grid that includes P ex t as a fourth parameter. 
In order to carry out option iii), the CMC model grid was 
produced for P ex t/kB = 10 4 cm -3 K (default), 10 5 ' 5 cm -3 K 
and 10 6 ' 5 cm -3 K, allowing for a look-up table of CO SLEDs 
in (Gq, togmc, Z', P ex t) parameter space. 

We examine the effects that options i)-iii) have on the 
global CO SLED of galaxy G2 (Figure 11) and how combi¬ 
nations of option i) and ii) change the values of the global 
CO-to-H 2 conversion factor for all three simulated galaxies 
(Table 4). Also, we show the impact that changes ii) and iii) 
have on the CO SLEDs of the individual CMC grid models 
in Figures C4 and C5. 

Changing the CMC mass spectrum from /3 = 1.8 to 
1.5 leaves the CO SLED virtually unchanged with only a 
marginal increase in the line flux ratios for J up > 4. This 
holds true regardless of what has been assumed for ii) and iii) 
(compare solid vs. dashed and dot-dashed vs. dotted curves 
in Figure 11). Also, changing /3 from 1.8 to 1.5 does not lead 
to any significant change in aco (Table 4). 

Adopting the modified Plummer density profile with 
an exponent of —7/2 instead of —5/2 results in significantly 
higher global line ratios for J up > 3 (see dotted vs. dashed 


Table 4. The global H 2 -to-CO conversion factors (in units of 
M@ pc -2 (Kkms -1 ) 1 ), averaged over Gl, G2, and G3, for com¬ 
binations of assumptions i) and ii). The pressure has been kept 
fixed at P ex t/^n = 10 4 cm -3 K. 



/3 = 1.8 

J = 1.5 

Plummer profile 

1.5 ±0.1 

1.5 ±0.1 

Modified Plummer profile 

3.6 ±0.4 

3.5 ±0.5 


curve and dot-dashed vs. solid curves in Figure 11). This is 
due to the higher central densities achieved for the modi¬ 
fied Plummer profile which, as shown in Figure C4, leads to 
an increase in the excitation of the higher J lines relative 
to that of the low-J lines. Significantly higher aco values 
(~ 3.6M 0 pc -2 (Kkms -1 ) -1 ) are obtained when adopting 
a modified Plummer profile (Table 4). These are in excel¬ 
lent agreement with the aco values inferred for 2 ~ 1.5 — 2 
BzK galaxies by Daddi et al. (2010). In our simulations, the 
higher aco values are due the steeper density profile which 
leads to smaller CMC sizes and therefore smaller CO(l — 0) 
fluxes (Figure C4) and, ultimately, a lowering of the total 
CO(l — 0) luminosity of the galaxy. 

Including external pressures of Pext/^B = 10 5 ' 5 cm 3 K 
and 10 6 ' 5 cm -3 K in the GMC model grid (see Section 3.4.4), 
means that these pressure values are assigned to GMCs with 
Pext/fe in the ranges 10 4 ' 75 —10 6 and > 10 e cm -3 K, respec¬ 
tively. The higher pressures now experienced by these sub¬ 
sets of GMCs results in smaller radii since -Rgmc oc P e x t -1 ^ 4 
for a given mass (see eq. 9). Their velocity dispersions and 
densities increase since a v oc P^ and rm 2 (R = 0) oc Pf/ 4 
for a given mass (see eqs. 10 and 11). The net effect this 
has on the global CO SLED is depicted in the bottom panel 
of Figure 11. Line ratios with J up = 3 are lowered, and in¬ 
creasingly so for higher J, resulting in a steeper decline at 
high (J up > 4) transitions. 


6 COMPARISON WITH OTHER MODELS 

The sub-grid physics implemented by SIGAME assumes that 
the scaling-laws and thermal balance equations that have 
been established for GMCs in our own Galaxy can be ap¬ 
plied to the ISM conditions in high -2 galaxies. In this re¬ 
spect, SIGAME does not differ from most other numerical 
simulations of the molecular line emission from galaxies 
(e.g., Narayanan et al. 2006; Pelupessy et al. 2006; Greve & 
Sommer-Larsen 2008; Christensen et al. 2012; Lagos et al. 
2012; Munoz & Furlanetto 2013; Narayanan & Krumholz 
2014; Lagos et al. 2012; Popping et al. 2014), although there 
are of course differences in the range and level of detail of 
the various sub-grid physics implementations. Here we high¬ 
light and discuss some of these differences between SIGAME 
and other simulations. 

First, however, we compare the sigame CO SLEDs of 
Gl, G2, and G3 with those predicted by other CO emis¬ 
sion simulations of similar main-sequence galaxies. A di¬ 
rect comparison can be made with the models presented by 
Narayanan & Krumholz (2014), where the global CO SLED 
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Figure 11. Global CO SLEDs of our model galaxy G2 for dif¬ 
ferent choices of ISM prescriptions. Top: For a pressure fixed of 
10 4 cm -3 K. Bottom: Using the P ex t as a fourth parameter in the 
GMC model grid (see Section 3.4.4). 


is parametrised as a function of the luminosity weighted 
SFR surface density (Esfr). We approximate a luminos¬ 
ity weighted Esfr by calculating a mean of Esfr for 1 kpc 2 
pixels across the galaxies seen face-on, weighted by the total 
SFR within the pixels. With this procedure, Gl, G2 and G3 
have Esfr of 3.0, 2.9 and 3.8 Mq yr _1 kpc 2 , respectively, 
and the resulting CO SLEDs inferred from the Narayanan 
& Krumholz (2014) parametrisation are shown as dotted 
lines in the bottom panel of Figure 10. These are seen to 
peak at J up = 5 and not 3 as the sigame CO SLEDs do. 
Also, the line flux ratios are significantly higher for all tran¬ 
sitions above J up = 2. Both set of models agree on the r 2 i 
ratio, but follow distinct CO SLED shapes for the remaining 
transitions. The models also seem to suggest one dominat¬ 
ing ISM phase, whereas the observed CO SLED of the BzK 
galaxies indicate a more complex multi-phased ISM. 

Combining a semi-analytical galaxy formation model 
with a photon-dominated region code, Lagos et al. (2012) 
made predictions of the CO SLEDs oi z = 2 galaxies as a 
function of their IR luminosities. They found that for galax¬ 
ies with infrared luminosities Lir ~ 10 11 ' 6 — 10 12:2 Lq - 
which is the range found for Gl, G2, and G3 if one adopts 
Lir/Lq = SFR/[Mq yr _1 ] x 10 1() (Kennicutt 1998 for a 
Chabrier IMF) - the CO SLEDs (when converted to ve¬ 


locity integrated flux ratios in order to compare with the 
bottom panel of Figure 10) peak at J = 3 — 2, in agree¬ 
ment with our model galaxies. However, the line ratios of 
the CO SLEDs of Lagos et al. (2012) are consistently above 
ours up to J up = 7, but cross over our models at J up = 7 
and go below at J up = 8. In terms of CO line luminosity, 
our model Gl agrees best with the models of Lagos et al. 
(2012), lying within the 10-90 percentile ranges (80% of their 
galaxy distribution) for the CO luminosities of their sample 
for 2 < J up < 8, but slightly below for the first two tran¬ 
sitions and for J up > 7. G2 and G3, lying roughly a factor 
2-4 above Gl, are outside the 10-90 percentile ranges of the 
models from Lagos et al. (2012). 

Popping et al. (2014), in their semi-analytical study of 
MS galaxies at z = 0, 1.2 and 2, found that for galaxies with 
far-infrared (FIR) luminosities of Lfir = 10 11 —10 12 Lq, the 
CO SLED peaks at CO(3—2) or CO(4—3) at z = 0, but at 
CO(6—5) at z = 2, as a result of denser and warmer gas in 
their model galaxies. The CO luminosities of our galaxies 
G2 and G3 are in broad agreement with the corresponding 
models by Popping et al. (2014) at CO(4—3) and (5—4) for 
similar SFRs. But at J up < 4 our galaxies are above while 
at J up > 5, our galaxies are below those of Popping et al. 
(2014) up until J up = 9 where the SLEDs cross again. 

• Implementation of G' 0 and £cr 

SIGAME stands out from most other simulations to date 
in the way the FUV radiation field and the CR flux that 
impinge on the molecular clouds are modelled and imple¬ 
mented. Most simulations adopt a fixed, galaxy-wide value 
of G' 0 and Ccr, scaled by the total SFR. or average gas 
surface density across the galaxy (e.g., Lagos et al. 2012; 
Narayanan & Krumholz 2014). SIGAME refines this scheme 
by determining a spatially varying G' 0 (and Ccr) set by the 
local SFRD, as described in Section 3.3. By doing so, we 
ensure that the molecular gas in our simulations is calori- 
metrically coupled to the star formation in their vicinity. 
If we adopt the method of Narayanan & Krumholz (2014) 
and calculate a global value for Go by calibrating to the MW 
value, using SFRmw = 2 Mq yr _1 and Go,mw = 0.6 Habing, 
our galaxies would have Go ranging from about 12 to 42, 
whereas, with the SFR surface density scaling of Lagos et al. 
(2012), global G' 0 values would lie between about 6 and 9. 
For comparison, the locally determined Go in our model 
galaxies spans a larger range from 0.3 to 27 (see Figure 
Cl). Narayanan & Krumholz (2014) determine the global 
value of £cr as 2 x 10 —17 ^ , s —1 , corresponding to values 
of 3.7 — 6.2 x 10~ 17 s _1 in our model galaxies, when us¬ 
ing the mass-weighted mean of Z' in each galaxy. Typi¬ 
cal values adopted in studies of ISM conditions are around 
(1 — 2) x 10 _17 s _1 (Wolfire et al. 2010; Glover & Clark 2012), 
but again, the local values of £cr in our galaxies span a 
larger range, from 1.3 x 10 -17 to 1.3 x 10 _15 s _1 . Adopt¬ 
ing the method of Narayanan & Krumholz (2014) for G' 0 
and £cr in our galaxies, leads to higher CO luminosities at 
J up > 2, but very similar CO(l — 0) luminosities which to¬ 
gether with slightly reduced molecular gas masses result in 
smaller aco factors by about 7 — 8 %. 

• The molecular gas mass fraction 

In SIGAME the molecular gas mass fraction (/(, lol ) is calcu¬ 
lated following the work by Pelupessy et al. (2006) (P06), in 
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Figure 12. Comparison of two methods for calculating the molec¬ 
ular gas mass fraction, f' j, of SPH particles in Gl: That of 
Pelupessy et al. (2006) with external cloud pressure derived from 
hydrostatic mid-plane equilibrium as done in this work (abscissa) 
and that of Krumholz et al. (2009) (ordinate) used by e.g., 
Narayanan & Krumholz (2014), color-coded by metallicity. 

which f' moX depends on temperature, metallicity, local FUV 
held, and boundary pressure on the gas cloud in question. 
Other methods exist, such as those of Blitz & Rosolowsky 
(2006) and Krumholz et al. (2009) (K09). The K09 method, 
which was adopted by e.g. Narayanan & Krumholz (2014), 
establishes /(n 0l from local cloud metallicity, dissociating ra¬ 
diation held and column density. In Figure 12 we compare 
the K09 method to that of P06. The gas surface density 
that enters in the K09 method was estimated within 1 kpc 
of each SPH particle, as used in our calculation of P ex t (see 
Section 3.3). 

The two methods agree for molecular gas fractions 
^ 80%, which are seen to be dominated by SPH particles 
with high metallicities (i.e., log Z' > 0.3). At lower molecu¬ 
lar gas fractions the agreement between the two methods ex¬ 
hibit an increasing scatter. Also, systematic, metallicity de¬ 
pendent differences are seen between the two methods. For 
high metallicity gas, the P06 method gives systematically 
higher molecular gas fractions than the K09 method. At low 
metallicities (i.e., \o%Z' £ 0.2) the K09 method tends to re¬ 
sult in higher molecular gas fractions than P06. The above is 
the result of the K09 method having a much weaker depen- 
dance on Z' Krumholz et al. (2009) than the P06 method 
(see Section 3.3). Despite these differences, we find that us¬ 
ing the K09 instead of the P06 method does not change 
the total molecular gas mass significantly: / mo i for Gl, G2 
and G3 goes from 12.2, 10.0, 9.2% to 11.4, 10.7 and 10.3%, 
respectively. 

A common feature of the above Hi —> H 2 prescriptions 
is that they assume an instantaneous Hi to H 2 conversion. 
However, as pointed out by Pelupessy et al. (2006), the typ¬ 
ical timescales for H 2 cloud formation are of order ~ 10' yr, 
which are comparable to the time-scales of a number of pro¬ 


cesses, such as cloud-cloud collisions and star formation, that 
can potentially alter or even fully disrupt typical molecular 
clouds and drive H 2 formation/destruction in both direc¬ 
tions. Ideally, therefore, a comprehensive modeling of the 
Hi —> H 2 transition should be time-dependent. Relying 
on the findings of Narayanan et al. (2011) and Krumholz 
& Gnedin (2011), the static solution for H 2 formation is a 
valid approximation for Z' > 0.01, which is the case for the 
three model galaxies studied here. 

• CO abundance 

SIGAME assumes a constant CO abundance relative to H 2 , 
which for the simulations presented in this paper was set to 
the Galactic CO abundance, i.e., [CO/H 2 ] = 2 x 10 -4 . In 
reality, H 2 gas can self-shield better than the CO gas, cre¬ 
ating an outer region or envelope of ‘CO-dark’ gas (Bolatto 
et al. 2008). Observations in the MW by Pineda et al. (2013) 
indicate that the amount of dark gas grows with decreasing 
metallicity and the resulting absence of CO molecules. By 
modeling a dynamically evolving ISM with cooling physics 
and chemistry incorporated on small scales, Smith et al. 
(2014) showed that in a typical MW-like disk of gas, the 
CO-dark gas mass fraction, /dg =, defined as having in¬ 
tegrated CO intensity Wco < O.IKkms -1 , is about 42%, 
but up to 62 % in a radiation field ten times that of the 
solar neighbourhood. While Smith et al. (2014) kept metal¬ 
licity fixed to solar, Wolfire et al. (2010) found /dg values of 
about 0.5 — 0.7 for the low metallicity cases with Z' = 0.5 in 
10 6 M@ clouds immersed in the local interstellar radiation 
field. Recently, Bisbas et al. (2015) argued that CR ioniza¬ 
tion rates of 10 — 50 x £cr,mw can effectively destroy most 
CO in GMCs. 

Thus, it is expected that [CO/H 2 ] is lower in regions 
of low metallicity and/or intense FUV and CR fields, ef¬ 
fectively leading to underestimates of the aco conversion 
factor when not accounted for in models. Since in this paper 
we have restricted our simulations to main-sequence galaxies 
with solar or higher than solar metallicities, adopting a con¬ 
stant Galactic [CO/H 2 ] seems a reasonable choice, but see 
Lagos et al. (2012), Narayanan et al. (2012) and Narayanan 
& Krumholz (2014) for alternative approaches. 

• CMC density and thermal structure 

When solving for the temperature structure of each CMC, 
SIGAME includes only the most dominant atomic and molec¬ 
ular species in terms of heating and cooling efficiencies. 
SIGAME also ignores heating via X-ray irradiation and tur¬ 
bulence. In particular, it has been shown with models that 
turbulent heating can completely control the mean tempera¬ 
ture of molecular clouds, exceeding the heating from cosmic 
rays (Pan & Padoan 2009). Our comparison of observed CO 
SLEDs at z = 2 with those of SIGAME, suggest that our sim¬ 
ulations are missing a warm gas component, which could 
be explained, at least in part, by the exclusion of turbulent 
heating. We have, however, checked our T k —uh 2 curves in 
Figure C2 against those of Glover & Clark (2012), who per¬ 
formed time-resolved, high-resolution (5m ~ 0.05 — 0.5 Mg) 
SPH simulations of individual FUV irradiated molecular 
clouds using a chemical network of 32 species (see their Fig¬ 
ure 2). Overall, there is great similarity in the Tk vs. riH 2 
behaviour of the two sets of simulations. This includes the 
same main trend of decreasing temperature with increas- 
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ing hydrogen density, as well as the local increase in Tk at 
~ 10 3 — 10 4 cm -3 and the subsequent decrease to Tk < 10 K 
at wh 2 > 10 5,5 cm -3 . Thus, our GMCs models, despite their 
simplified density profiles and chemistry, seem to agree well 
with much more detailed simulations. 


7 SUMMARY 

In this paper we have presented SIGAME, a code that simu¬ 
lates the molecular line emission of galaxies via a detailed 
post-processing of the outputs from cosmological SPH sim¬ 
ulations. A sequence of sub-grid prescriptions are applied 
to a simulation snapshot in order to derive the molecular 
gas density and temperature from the SPH particle infor¬ 
mation, which includes SFR, gas density, temperature and 
metallicity. 

A key aspect of SIGAME is the localised coupling between 
the star formation and the energetics of the ISM, where the 
strength of the local FUV radiation field and CR ionization 
rate that impinge and heat a cloud scales with the local 
star formation rate density. The radial temperature profile 
of each GMC is calculated by balancing the heating rate 
with cooling from H 2 , CO, [Ol], and [Cii] lines in addition 
to gas-dust interactions. In these thermal balance calcula¬ 
tions, SIGAME takes into account the local enrichment of the 
gas, which is a critical parameter for the cooling rates. The 
CO emission line spectrum from a grid of GMC models is 
calculated using the 3D radiative transfer code lime. 

We used SIGAME to create line emission velocity-cubes 
of the full CO rotational ladder for three cosmological N- 
body/SPH simulations of massive (M* > 10 10,5 Mg) main- 
sequence galaxies at 2 = 2. 

Molecular gas is produced more efficiently towards the 
centre of each galaxy, and while Hi surface gas densities (in¬ 
cluding helium) do not exceed ~ 100 M@ pc -2 anywhere 
in the disk, central molecular gas surface densities reach 
~ 1000 Mq pc -2 on spatial scales of 100 pc x 100 pc, in 
good agreement with observations made at similar spatial 
resolution. This strong increase in molecular surface den¬ 
sity is brought on by a similar increase in total gas surface 
density, overcoming the increase in photo-dissociating FUV 
held towards the centre of each galaxy. 

Turning to the CO emission, the velocity-integrated mo¬ 
ment 0 maps reveal distinct differences in the various transi¬ 
tions as molecular gas tracers. The morphology of molecular 
gas in our model galaxies is well reproduced in CO (1-0), but 
going to higher transitions, the region of CO emitting gas 
shrinks towards the galaxy centres. The global CO SLEDs of 
our simulated galaxies all peak at J = 3—2. Recent CO(5—4) 
observations of 2 ~ 1.5 BzK galaxies seem to suggest that 
these galaxies actually peak at higher J. The CO(3 —2) line 
luminosities of our model galaxies are within the range of 
corresponding observed samples at redshifts 2 ~ 1 — 2.5, 
however on the low side. In particular, the model galaxies 
are below or at the CO luminosities of BzK-selected galaxies 
of comparable mass and SFR but at 2 ~ 1.5. The low lumi¬ 
nosities are most likely a consequence of molecular gas mass 
fractions in our galaxies being about ~ 4 — 5 times below the 
observed values in the star-forming galaxies at 2 = 1 — 2.5 
used to compare with. 


Combining the derived H 2 gas masses with the CO(l — 
0) line emission found, we investigate local variations in the 
CO-H 2 conversion factor aco- The radial otco profiles all 
show a decrease towards the galaxy centres, dropping by a 
factor of ~ 1.2 in the central R < 2kpc region compared to 
the disk average, the main driver being the FUV field rather 
than a gradient in density or metallicity. Global aco factors 
range from 1.4 to 1.6Mg pc -2 (Kkms -1 ) -1 or about 0.3 
times the MW value, but closer to values for 2 ~ 1.5 normal 
star-forming galaxies identified with the BzK colour criteria. 
Changing the GMC properties from what is observed in the 
MW and local galaxies to a steeper GMC density profiles 
and/or a shallower GMC mass spectra, resulted in elevated 
aco values of up to 0.9 times the MW value. 

The CO luminosity ratios of CO 3—2/1—0 and 
7—6/1—0 (r 32 and r^e respectively) drop off in radius about 
where the FUV radiation drops in intensity, and are thus 
likely controlled by the FUV field as is aco- The global 
ratios of r 2 i — 1 and r 32 — 0.6 agree very well with observa¬ 
tions of BzK galaxies, while the rs 4 of about 0.15 is low com¬ 
pared to recent observations in BzK—4171 and BzK—21000. 
The low flux from high CO transitions in our models com¬ 
pared to observations could be explained, at least in part, by 
our omission of turbulent heating in SIGAME. However, more 
observations of J up > 3 lines towards high -2 main-sequence 
galaxies, such as the BzKs, are still needed in order to deter¬ 
mine the turn-over in their CO SLEDs and better constrain 
the gas excitation. 

Finally, we note that SIGAME in principle is able to sim¬ 
ulate the emission from a broad range of molecular and 
atomic lines in the far-IR/mm wavelength regime provided 
measured collision rates exist, such as those found in the 
LAMDA database 4 . In the future, we aim to apply sigame 
to a larger and more diverse sample of galaxies, possibly at 
higher redsliift, and explore the emission from other impor¬ 
tant ISM diagnostic lines. 


ACKNOWLEDGMENTS 

We are grateful to Inti Pelupessy for elaborations on his pa¬ 
pers and help with the ionization fractions, and to Desika 
Narayanan for stimulating discussions. We thank Padelis 
Papadopoulos for useful discussions and suggestions along 
the way. We also thank Julia Kamenetzky for providing us 
with CO line data from her paper. The draft was much im¬ 
proved by thoughful reports from the anonymous referee. 
KPO gratefully acknowledge the support from the Lundbeck 
foundation and TRG acknowledges support from a STFC 
Advanced Fellowship. S. T. acknowledges support from the 
ERC Consolidator Grant funding scheme (project ConTExt, 
grant number 648179). The Dark Cosmology Centre and 
Starplan are funded by the Danish National Research Foun¬ 
dation. 

REFERENCES 

Ackermann M., et al., 2013, Science, 339, 807 

4 http://www.strw.leidenuniv.nl/-moldata/, Schoier et al. 
(2005) 


MNRAS 000, 1-?? (2015) 


Simulating CO emission with SIGAME 21 


Adelberger K. L., Steidel C. C., Shapley A. E., Hunt M. P., Erb 
D. K., Reddy N. A., Pettini M., 2004, ApJ, 607, 226 
Aravena M., et al., 2010, ApJ, 718, 177 
Aravena M., et al., 2014, MNRAS, 442, 558 
Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 
47, 481 

Bakes E. L. O., Tielens A. G. G. M., 1994, ApJ, 427, 822 
Barvainis R., Tacconi L., Antonucci R., Alloin D., Coleman P., 
1994, Nature, 371, 586 

Bigiel F., Leroy A., Walter F., Brinks E., de Blok W. J. G., 
Madore B., Thornley M. D., 2008, AJ, 136, 2846 
Black J. H., Dalgarno A., 1977, ApJs, 34, 405 
Blanc G. A., et al., 2013, ApJ, 764, 117 
Blitz L., Rosolowsky E., 2006, ApJ, 650, 933 
Blitz L., Fukui Y., Kawamura A., Leroy A., Mizuno N., 
Rosolowsky E., 2007, Protostars and Planets V, pp 81-96 
Bolatto A. D., Leroy A. K., Rosolowsky E., Walter F., Blitz L., 
2008, ApJ, 686, 948 

Bolatto A. D., Wolfire M., Leroy A. K., 2013, ARA&A, 51, 207 
Bothwell M. S., et al., 2013, ApJ, 779, 67 
Bovy J., Rix H.-W., Hogg D. W., 2012, ApJ, 751, 131 
Brinch C., Hogerheijde M. R., 2010, A&A, 523, A25 
Bryant P. M., Scoville N. Z., 1999, AJ, 117, 2632 
Carilli C. L., Walter F., 2013, ARAA, 51, 105 
Carroll T. J., Goldsmith P. F., 1981, ApJ, 245, 891 
Cazaux S., Spaans M., 2004, ApJ, 611, 40 
Chabrier G., 2003, PASP, 115, 763 
Chomiuk L., Povich M. S., 2011, AJ, 142, 197 
Christensen C., Quinn T., Governato F., Stilp A., Shen S., Wad- 
sley J., 2012, MNRAS, 425, 3058 
Daddi E., Cimatti A., Renzini A., Fontana A., Mignoli M., 
Pozzetti L., Tozzi P., Zamorani G., 2004, ApJ, 617, 746 
Daddi E., et al., 2010, ApJ, 713, 686 
Daddi E., et al., 2015, A&A, 577, A46 

Dame T. M., Hartmann D., Thaddeus P., 2001, ApJ, 547, 792 
Dannerbauer H., Daddi E., Riechers D. A., Walter F., Carilli 
C. L., Dickinson M., Elbaz D., Morrison G. E., 2009, ApJL, 
698, L178 

Downes D., Solomon P. M., 1998, ApJ, 507, 615 
Draine B. T., 2011, Physics of the Interstellar and Intergalactic 
Medium. Princeton University Press 
Elmegreen B. G., 1989, ApJ, 344, 306 

Fixsen D. J., Bennett C. L., Mather J. C., 1999, ApJ, 526, 207 
Frayer D. T., Ivison R. J., Scoville N. Z., Yun M., Evans A. S., 
Smail I., Blain A. W., Kneib J.-P., 1998, ApJL, 506, L7 
Geach J. E., Papadopoulos P. P., 2012, ApJ, 757, 156 
Glover S. C. O., Clark P. C., 2012, MNRAS, 421, 9 
Gnat O., Sternberg A., 2007, ApJS, 168, 213 
Goldsmith P. F., 2001, ApJ, 557, 736 
Greve T. R., Sommer-Larsen J., 2008, A&A, 480, 335 
Haardt F., Madau P., 2001, preprint, (arXiv:astro-ph/0106018) 
Heiderman A., Evans II N. J., Allen L. E., Huard T., Heyer M., 
2010, ApJ, 723, 1019 

Heyer M. H., Brunt C. M., 2004, ApJL, 615, L45 
Heyer M., Dame T. M., 2015, ARA&A, 53, 583 
Hughes A., et al., 2013, ApJ, 779, 44 
Iono D., et al., 2009, ApJ, 695, 1537 
Jarrett T. H., et al., 2013, AJ, 145, 6 

Kamenetzky J., Rangwala N., Glenn J., Maloney P. R., Conley 
A., 2015, preprint, (arXiv: 1508.05102) 

Kennicutt Jr. R. C., 1998, ARA&A, 36, 189 
Krumholz M. R., Gnedin N. Y., 2011, ApJ, 729, 36 
Krumholz M. R., McKee C. F., Tumlinson J., 2009, ApJ, 699, 
850 

Lagos C. d. P., Bayet E., Baugh C. M., Lacey C. G., Bell T. A., 
Fanidakis N., Geach J. E., 2012, MNRAS, 426, 2142 
Larson R. B., 1981, MNRAS, 194, 809 

Lee H.-H., Bettens R. P. A., Herbst E., 1996, A&As, 119, 111 


Leroy A. K., Walter F., Brinks E., Bigiel F., de Blok W. J. G., 
Madore B., Thornley M. D., 2008, AJ, 136, 2782 
Leroy A. K., et al., 2015, ApJ, 801, 25 
Magnelli B., et al., 2012, A&A, 548, A22 

Mao R.-Q., Schulz A., Henkel C., Mauersberger R., Muders D., 
Dinh-V-Trung 2010, ApJ, 724, 1336 
McMillan P. J., 2011, MNRAS, 414, 2446 
Mentuch Cooper E., et al., 2012, ApJ, 755, 165 
Monaghan J. J., 2005, Reports on Progress in Physics, 68, 1703 
Munoz J. A., Furlanetto S. R., 2013, MNRAS, 435, 2676 
Narayanan D., Hopkins P. F., 2013, MNRAS, 433, 1223 
Narayanan D., Krumholz M. R., 2014, MNRAS, 442, 1411 
Narayanan D., et al., 2006, ApJL, 642, L107 
Narayanan D., et al., 2008a, ApJS, 174, 13 
Narayanan D., et al., 2008b, ApJs, 176, 331 

Narayanan D., Cox T. J., Shirley Y., Dave R., Hernquist L., 
Walker C. K., 2008c, ApJ, 684, 996 
Narayanan D., Cox T. J., Hayward C. C., Younger J. D., Hern¬ 
quist L., 2009, MNRAS, 400, 1919 
Narayanan D., Krumholz M., Ostriker E. C., Hernquist L., 2011, 
MNRAS, 418, 664 

Narayanan D., Krumholz M. R., Ostriker E. C., Hernquist L., 
2012, MNRAS, 421, 3127 
Neri R., et al., 2003, ApJL, 597, LI 13 
Norman C. A., Spaans M., 1997, ApJ, 480, 145 
Obreschkow D., Klockner H.-R., Heywood I., Levrier F., Rawlings 
S., 2009, Apj, 703, 1890 

Obreschkow D., Heywood I., Rawlings S., 2011, Apj, 743, 84 
Pan L., Padoan P., 2009, ApJ, 692, 594 

Papadopoulos P. P., Thi W.-F., 2013, in Torres D. F., Reimer O., 
eds, Advances in Solid State Physics Vol. 34, Cosmic Rays in 
Star-Forming Environments. Springer-Verlag, Berlin, Heidel¬ 
berg, p. 41 (arXiv: 1207.2048), doi: 10.1007/978-3-642-35410- 
6'5 

Papadopoulos P. P., Isaak K., van der Werf P., 2010, Apj, 711, 
757 

Papadopoulos P. P., Thi W.-F., Miniati F., Viti S., 2011, MNRAS, 
414, 1705 

Papadopoulos P. P., van der Werf P. P., Xilouris E. M., Isaak 
K. G., Gao Y., Muhle S., 2012, MNRAS, 426, 2601 
Papadopoulos P. P., et al., 2014, ApJ, 788, 153 
Papovich C., Finkelstein S. L., Ferguson H. C., Lotz J. M., Gi- 
avalisco M., 2011, MNRAS, 412, 1123 
Pelupessy F. I., 2005, PhD thesis, Leiden Observatory, Leiden 
University, P.O. Box 9513, 2300 RA Leiden, The Netherlands 
Pelupessy F. I., Papadopoulos P. P., van der Werf P., 2006, ApJ, 
645, 1024 

Pineda J. E., Caselli P., Goodman A. A., 2008, Apj, 679, 481 
Pineda J. L., Langer W. D., Velusamy T., Goldsmith P. F., 2013, 
A&A, 554, A103 

Planck Collaboration XVI 2014, A&A, 571, A16 
Plummer H. C., 1911, MNRAS, 71, 460 

Popping G., Perez-Beaupuits J. P., Spaans M., Trager S. C., 
Somerville R. S., 2014, MNRAS, 444, 1301 
Price D. J., 2007, Publications of the Astronomical Society of 
Australia, 24, 159 

Puchwein E., Bolton J. S., Haehnelt M. G., Madau P., Becker 
G. D., Haardt F., 2015, MNRAS, 450, 4081 
Riechers D. A., et al., 2011, ApJL, 739, L32 
Rodighiero G., et al., 2010, A&A, 518, L25 

Rollig M., Ossenkopf V., Jeyakumar S., Stutzki J., Sternberg A., 
2006, A&A, 451, 917 

Romeo A. D., Sommer-Larsen J., Portinari L., Antonuccio-Delogu 
V., 2006, MNRAS, 371, 548 
Saintonge A., et al., 2011, MNRAS, 415, 32 
Sandstrom K. M., et al., 2013, ApJ, 777, 5 
Schoier F. L., van der Tak F. F. S., van Dishoeck E. F., Black 
J. H., 2005, A&A, 432, 369 


MNRAS 000, 1-?? (2015) 


22 K. P. Olsen 


Schuster K. F., Kramer C., Hitschfeld M., Garcia-Burillo S., 
Mookerjea B., 2007, A&A, 461, 143 
Seon K.-L, et al., 2011, ApJs, 196, 15 
Simnett G. M., McDonald F. B., 1969, ApJ, 157, 1435 
Smith B., Sigurdsson S., Abel T., 2008, MNRAS, 385, 1443 
Smith R. J., Glover S. C. O., Clark P. C., Klessen R. S., Springel 
V., 2014, MNRAS, 441, 1628 

Sofia U. J., Lauroesch J. T., Meyer D. M., Cartledge S. I. B., 
2004, ApJ, 605, 272 

Solomon P. M., Downes D., Radford S. J. E., Barrett J. W., 1997, 
ApJ, 478, 144 

Sommer-Larsen J., Vedel H., Hellsten U., 1998, MNRAS, 294, 485 
Sommer-Larsen J., Gotz M., Portinari L., 2003, ApJ, 596, 47 
Sommer-Larsen J., Romeo A. D., Portinari L., 2005, MNRAS, 
357, 478 

Speagle J. S., Steinhardt C. L., Capak P. L., Silverman J. D., 
2014, ApJS, 214, 15 

Springel V., Hernquist L., 2003, MNRAS, 339, 289 
Stahler S. W., Palla F., 2005, The Formation of Stars. Wiley 
Stinson G., Seth A., Katz N., Wadsley J., Governato F., Quinn 
T., 2006, MNRAS, 373, 1074 
Strong A. W., Mattox J. R., 1996, A&A, 308, L21 
Swinbank A. M., et al., 2011, ApJ, 742, 11 
Tacconi L. J., et al., 2008, ApJ, 680, 246 
Tacconi L. J., et al., 2013, ApJ, 768, 74 

Tielens A. G. G. M., 2005, The Physics and Chemistry of the 
Interstellar Medium. Cambridge University Press 
Vlahakis C., van der Werf P., Israel F. P., Tilanus R. P. J., 2013, 
MNRAS, 433, 1837 
Webber W. R., 1998, ApJ, 506, 329 

Weifi A., Downes D., Neri R., Walter F., Henkel C., Wilner D. J., 
Wagg J., Wiklind T., 2007, A&A, 467, 955 
Whitaker K. E., et al., 2011, ApJ, 735, 86 
Whitworth A. P., Ward-Thompson D., 2001, ApJ, 547, 317 
Wiersma R. P. C., Schaye J., Smith B. D., 2009, MNRAS, 393, 
99 

Wolfire M. G., McKee C. F., Hollenbach D., Tielens A. G. G. M., 
2003, ApJ, 587, 278 

Wolfire M. G., Hollenbach D., McKee C. F., 2010, ApJ, 716, 1191 
Wright E. L., et al., 1991, ApJ, 381, 200 
Wu R., et al., 2015, A&A, 575, A88 

Yang B., Stancil P. C., Balakrishnan N., Forrey R. C., 2010, ApJ, 
718, 1062 

Yun M. S., et al., 2015, MNRAS, 454, 3485 
Zahid H. J., Dima G. I., Kewley L. J., Erb D. K., Dave R., 2012, 
ApJ, 757, 54 


APPENDIX A: THERMAL BALANCE OF THE 
ATOMIC GAS PHASE 

As explained in § 3.2, we cool the initial hot SPH gas by 
requiring: Tcr ,Hl — Aatoms+ions + A 

rec + Af_f. 

Tcr.hi is the heating rate of the atomic gas due to cos¬ 
mic ray ionizations (Draine 2011): 


rcR,m =1.03 x 10 "'ram 


1 + 4.06 


( CCR.HI 

VlfFie 

Xe 

X e + 0.07 


1 / 2 ' 


-3 -1 
erg cm s , 


(Al) 


where Ccr.hi is the primary CR ionization rate of Hi atoms 
(determined locally in our simulations according to eq. 12), 
and x e is the hydrogen ionization fraction calculated with 


a procedure kindly provided by I. Pelupessy; see also Pelu- 
pessy (2005). The term containing x e in eq. Al accounts for 
the fact that in highly ionised gas, electrons created by pri¬ 
mary CR ionization have a high probability of transferring 
their kinetic energy into heat via long-range Coulomb scat¬ 
tering off free electrons. For low ionization gas, this term 
becomes insignificant as a higher fraction of the energy of 
the primary electrons goes to secondary ionizations or exci¬ 
tation of bound states instead of heating. 

A a toms+ions is the total cooling rate due to line emission 
from H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe, calcu¬ 
lated using the publically available code of Wiersma et al. 
(2009) which takes T k , riH and the abundances of the above 
elements as input. Wiersma et al. (2009) compute the cool¬ 
ing rates with the photoionization package CLOUDY assuming 
CIE. They also adopt a value for the meta-galactic UV and 
X-ray field equal to that expected at z ~ 2 (Haardt & Madau 
2001). At z ~ 2, the emission rate of Hi ionizing radiation 
is higher by a factor of about ~ 30 than at z = 0 (Puchwein 
et al. 2015), and thus plays an important role in metal line 
cooling calculations. 

A rec is the cooling rate due to hydrogen recombination 
emission (Draine 2011): 

A rec = aBn e n H + (E rr ) ergscm _3 s _1 , (A2) 

where as is the radiative recombination rate for hydrogen 
in the case of optically thick gas in which ionizing photons 
emitted during recombination are immediately re-absorbed. 
We adopt the approximation for as given by Draine (2011): 

ob = 2.54 X 10 T 4 ' 4; cm s , (A3) 

where T 4 is defined as Tk/10 4 K. The density of ionised hy¬ 
drogen, n H +, is set equal to the electron density, n e , and E rr 
is the corresponding mean kinetic energy of the recombining 
electrons: 

(Err) = [0.684 - 0.0416 lnP 4 ] k B T k ergs. (A4) 

Af_f is the cooling rate due to free-free emission from 
electrons in a pure H plasma (i.e., free electrons scattering 
off H + ), and is given by (Draine 2011): 

Af_f = 0.54TJ' 37 fcflTkRen H +aB ergscm _3 s _1 , (A5) 

where the recombination rate, a B , is calculated in the same 
way as for A rec . 

Figure Al shows the above heating and cooling rates 
pertaining to two example SPH particles with similar initial 
temperatures (~ 10 4 K). Because of different ambient condi¬ 
tions (i.e., nm, x e , Z' , and Ccr) the equilibrium temperature 
solutions for the two gas particles end up being significantly 
different. 


APPENDIX B: THERMAL BALANCE OF THE 
MOLECULAR GAS PHASE 

As described in Section 3.4 SIGAME assumes that the molec¬ 
ular gas resides exclusively in giant molecular clouds that 
have Plummer radial density profiles (i.e., given by eq. 11). 
Throughout the clouds the gas temperature is solved for ac¬ 
cording to the heating and cooling equilibrium requirement 
r pe T r cr,h 2 = Ah 2 +Aco+Acii + Aoi-l-Ag a ,j_dust (eq. 14). 
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Figure Al. Heating and cooling rates as functions of temperature for two SPH gas particles with [npj = 14.62 cm 3 , = 4.02 x 

10 -16 s -1 , :r e (T k , S PH = 8892 K) = 0.006] (left) and [n H = 0.09 cm- 3 , Ccr = 3.92 x lO -1 ^ -1 , z e (T k , S PH = 5649 K) = 0.001] (right). 
In both cases the metal line emission is the main cooling agent. The left-hand side plot illustrates the case of an SPH particle of high gas 
density and metallicity, leading to relatively efficient metal line cooling and in a low equilibrium temperature of T k = 37 K. The right- 
hand side plot shows a case of lower density and metallicity causing less cooling by metal lines, hence a higher equilibrium temperature 
despite a slightly lower CR heating rate. The dashed vertical lines in the two panels indicate the original SPH gas temperatures of the 
SPH particles, and the solid vertical lines mark their final equilibrium temperatures. 


Tpe is the heating rate of the gas due to photo-electric 
ejection of electrons from dust grains by FUV photons, and 
is given by (Bakes &; Tielens 1994): 

r PE = 10" 24 eGo, att n H ergscm -3 s -1 , (Bl) 


where G^att is the local attenuated FUV held in Habing 
units, derived following eq. 13, and e is the heating efficiency: 


4.87 x 10“ 


l + 4x 10- 3 (Gi iatt T 0 - 5 /ne) 0 - 73 
3.65 x 10 -2 

+ l + 4x 10- 3 (G;,, att TO'Vne) 0 - 73 ’ 


(B2) 


where n e is the electron density, calculated as x e nn, with x e 
again calculated using the procedure of I. Pelupessy. 

Pcr.h, is the heating rate by cosmic rays traveling 
through molecular gas (Stahler & Palla 2005): 


Pcr,h 2 = 1.068xl0 -24 


( Ccr,h 2 \ 

v 10 -16 ) 


( n H2 \ 

V10 3 cm -3 / 


-3 -1 

ergs cm s 


(B3) 


where Ccr,h 2 is the local CR primary ionization rate of H 2 
molecules, which is approximately 1.6 x higher that of Hi 
atoms (Stahler & Palla 2005). 

Ah 2 is the H 2 line cooling rate, and we use the parame¬ 
terization made by Papadopoulos et al. (2014) that includes 
the two lowest H 2 rotational lines (S(0) and S(l), the only 
lines excited for Tk < 1000 K): 


Ah 2 =2.06 x 10 


— 24 ^H2 


1 + r Q 


1 4. I J510K/T k 

+ 5 


x (1 + Rio) ergs cm 3 s 1 , 



(B4) 


where Rio is defined as: 


Rio = 26.8 r op 


1 + (l/5)e 510K/Tk (l + 4 Hl) ' 
l + (3/7)e8«K/T k (l + ^)_ ’ 


(B5) 


and no ~ 54 cm -3 and ni ~ 10 3 cm -3 are the critical den¬ 
sities of the S(0):2-0 and S(l):3-1 rotational lines. r op is the 
ortho-H 2 /para-H 2 ratio (set to 3 which is the equilibrium 
value). 

For the cooling rates due to the [Cll]158/im and 
[0 1 ] 63 /im+146 /im fine-structure lines we adopt the param- 
eterizations by Rollig et al. (2006). The Cll cooling rate 
(Acn) is: 


Aon =2.02 x 10 - 24 nZ' 


(B6) 


1 + Ie 92K / Tk (1 + 1300/n H ) 


-3 -1 
ergs cm s , 


where a carbon to hydrogen abundance ratio that scales with 
metallicity according to X[C] = 1-4 x 10 -4 Z' is assumed. For 
the parameterization of the 0 1 cooling rate (Aoi = A 63 Mm + 
Ai 46 M m) we refer to eqs. A.5 and A.6 in Rollig et al. (2006) 
and simply note that we adopt (in accordance with Rollig 
et al. 2006) an oxygen to hydrogen abundance ratio of X[o] = 
3 x 10 ~ A Z'. 

Aco is the cooling rate due to CO rotational transitions. 
We use the parameterization provided by Papadopoulos & 
Thi (2013): 


a iivin- 24 ( nH! ) 3/ V V fxc°\ -3 -1 

A CO = 4.4x10 (—Jergscm s , 

(B7) 

where Xco/X[C] is the relative CO to neutral carbon abun¬ 
dance ratio, the value of which we determine by interpo¬ 
lation, assuming that Xco/X[C] = (0.97,0.98,0.99,1.0) for 
tih 2 = 5 x 10 3 ,10 4 ,10 5 ,10 6 cm -3 , respectively (Papadopou¬ 
los & Thi 2013). 

Ag a s_dust is the cooling rate due to gas-dust interactions 
and is given by (Papadopoulos et al. 2011): 

Agas-dust = 3.47 x 10 - 33 n H 2 \/21(Tk - T dust ) ergs cm -3 s -1 , 

(B8) 
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where the dust temperature (Tdust) is calculated using eq. 15 
(Section 3.4.3). 


APPENDIX C: GMC MODELS 

Each SPH particle is divided into several GMCs as described 
in § 3.4.1, and we derive the molecular gas density and tem¬ 
perature within each from three basic parameters which are 
SFR density, GMC mass, mGMC, and metallicity, Z' /Zq. 
Derived from these basic parameters are the far-UV and cos¬ 
mic ray field strengths, the Ha gas mass fraction of each SPH 
particles, as well as the GMC properties used to derive the 
CO excitation and emission; H 2 density and temperature. 
Histograms of the basic parameters are shown in Figured, 
while properties derived thereof can be found in Figure C2. 

This paper has been typeset from a TpX/PTpX file prepared by 
the author. 
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Figure Bl. Equilibrium heating (red curves) and cooling (blue curves) rates of the gas as functions of H 2 density for two different GMC 
models with [mcMC = 10 4 M©, Z' = 1, Gq = 27, £cr = 1.35 X 10 -15 s -1 , P ex t = 10 4 Kcm -3 ] (left) and [mcMC = 1O 4 M0, Z' = 63, 
Gq = 1, £cr = 5 x 10 -17 s -1 , Pext = 10 4 Kcm -3 ] (right). In the case of high FUV and CR radiation fields but low metallicity (left), 
the heating is dominated by cosmic ray heating (red solid) in the inner region (nn 2 > 1500cm -3 ) and photoelectric heating (red dotted) 
in the outer region. Cooling is dominated by gas-dust interactions (blue dashed) in the inner region ( nn 2 > 10000 cm -3 ) and by [Cn] as 
well as H 2 line cooling (blue long-dashed and dotted) in the outer region. In the opposite case of low FUV and CR radiation fields but 
high metallicity (right), the same heating and cooling mechanisms are dominating the energy balance throughout the cloud, except in 
the inner region, where cooling by CO line emission (blue solid) is more important than it is in the case of low metallicity. 
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Figure Cl. Mass-weighted histograms of the basic parameters of the GMCs in G1 (red solid), G2 (green dashed), and G3 (blue dotted). 
From top left and clockwise: the local far-UV field (Gq) (and CR ionization rate since (cr oc Gq), GMC mass (wgmc)i metallicity ( Z ') 
and external pressure (P e xt)- Black vertical lines indicate the Gq, mcMCi Z' 5 ^ext -values for which — nn 2 curves were calculated (see 
Figure C2) — a total of 630 GMCs which make up our grid GMC models. Each GMC in the galaxies is assigned the — nn 2 curve of 
the GMC model at the closest grid point. 
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Figure C2. Kinetic temperature versus H2 density curves for 80 out of the 630 grid model GMCs that span the full (Gq, Mqmc? Z ') 
parameter space set by and marked on top of the distributions in Figure Cl (see also Section 3.4.4) for a pressure of P e xt = 10 4 Kcm -3 . 
The grid model most often assigned to GMCs in G1 is indicated by the red dashed curve in the highlighted panel and corresponds to 
G' 0 = 7.0 (Ccr — 3.5 x 10 —16 s —1 ), logmQMc/M© = 4.25, and Z' = 3.2. In general, higher metallicity (from top to bottom) leads to 
more cooling via emission lines of ions, atoms and molecules, and hence lower temperatures. On the other hand, higher UV and CR 
fields (from left to right) cause more heating and therefore higher T^. The decreasing trend of with higher values of nn 2 is mainly 
caused by the gradual attenuation of the UV field as one moves into the cloud. The ‘bump’ at tih 2 10 3 — 10 4 cm -3 corresponds to 
the transition from CII line cooling to the less efficient CO line cooling at higher densities. At densities above nn 2 = 10 4 cm -3 , gas-dust 
interactions set in and eventually cools the gas down to the CMB temperature in all GMC cores. 
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Figure C3. CO SLEDs obtained with lime for the 80 model GMCs whose -nn 2 curves are shown in Figure C2. The CO SLED used 
most often in G1 is shown as the red dashed curve in the highlighted panel. These CO SLEDs were made for a fixed external pressure 
of -Pext/&B — 10 4 cm -3 K and a default Plummer density profile. 
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Figure C4. CO SLEDs obtained with lime for the same \G'q,Z'\ values as shown in the panels of Figure C3 for Plummer density 
profiles with power-law index —5/2 (solid curve) and —7/2 (dashed curve). In all panels, the external pressure has been fixed to 
Pext/^B — 10 4 cm -3 K and the GMC mass to mcMC — 10 4 ' 25 Mq. 
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Figure C5. CO SLEDs obtained with lime for the same [GqjZ'] values as shown in the panels of Figure C3 for Pext/^B = 10 4 cm -3 K 
(solid curves), 10 5 5 cm -3 K (dashed curves), and 10 6,5 cm -3 K (dot-dashed curves). In all panels, the GMC mass is fixed to 10 4 25 Mq. 
Higher pressure environments are seen to lead to a decrease in luminosity for all transitions. 
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APPENDIX D: TESTING sigame ON MW-LIKE 
GALAXIES 

As a benchmark test, SIGAME was applied to galaxy simulations 
at z = 0 with properties similar to those of the MW. The reason¬ 
ing here was that since much of the sub-grid physics adopted by 
SIGAME relies on empirical relations observed in the MW (e.g., the 
GMC mass spectrum index), the method ought to come close to 
re-producing the observed CO properties of spiral galaxies, such 
as our MW, in the local universe. 

The three galaxies (hereafter denoted MW1, MW2, and 
MW3 in order of increasing SFR) were selected from a more 
recent version of the SPH simulation described in Section 2. 
The properties of MW1, MW2, and MW3 are listed in Ta¬ 
ble D1 and are seen to be within a factor ~ 2 — 6 of the stel¬ 
lar mass (M* ? mw — 6 x 10 10 Mq; McMillan 2011) and SFR 
(1.9 zb 0.4 Mq yr -1 ; Chomiuk & Povich 2011) of the MW. 

The steps outlined in Section 3.1 (and further described 
in Sections 3.2 to 3.4.3) were followed, but with two minor 
modifications: 1) the CMB temperature at z = 0 was set to 
2.725 K rather than 8.175K used at z = 2, and 2) the resolved 
FUV fields of our MW-like model galaxies span a range below 
that seen in our z = 2 model galaxies (cf. Figured) meaning 
that the Gq parameter values for our GMC grid had to be ad¬ 
justed correspondingly. The following Gq grid points were chosen: 
[0.05, 0.1, 0.15, 0.1, 0.2, 0.4, 0.6, 0.8,1.0,1.2] Habing. 

In Figure D1 we compare the resulting CO SLEDs with that 
of the MW from Fixsen et al. (1999), as well as with those 
of M51 (M* = 4.7 x 10 lo M©; Mentuch Cooper et al. (2012), 
SFR = 2.6 M© yr —1 ; Schuster et al. 2007) from Vlahakis et al. 
(2013), Hughes et al. (2013) and Kamenetzky et al. (2015), M83 
(M* = 7.9 x 10 10 M© and SFR = 3.2 M© yr - - 1 ; Jarrett et al. 
2013) from Wu et al. (2015) and other local galaxies of IR lu¬ 
minosity within 20% of that of the MW (1.8 x 10 10 L©; Wright 
et al. 1991) from Kamenetzky et al. (2015). M51 is a nearby 
galaxy slightly smaller in size and stellar mass than the MW, 
whereas M83 is more massive and star-forming than the MW. 
The dispersion in observed CO SLEDs among the MW, M83, 
M51 and other local galaxies, allows for a large range in nor¬ 
malization and shape of the CO SLED within which our model 
galaxies find themselves. 

The CO line luminosities of the MW lie roughly in the range 
spanned by the CO SLEDs of MW2 and MW3 up to and including 
J up = 4, but significantly above that of MW1. For the higher J 
transitions the line luminosities of our simulations drop off more 
rapidly than the MW, meaning that our model galaxies have lower 
luminosities compared to the MW (and other local galaxies of IR 
luminosities within 20% from that of the MW; long-dashed orange 
lines) at J up > 4. This suggests that our simulations are missing 
a warm/dense component which is required to significantly excite 
these high-J lines. The agreement between the CO SLED of MW2 
and that of M51 is good at low transitions (J up < 6), which is 
encouraging given that the two galaxies have nearly identical SFR 
and M*. At higher J values, the CO SLED of M51 has an excess 
of CO emission when compared with MW2, corresponding to a 
warm gas phase not captured by our models. This warm phase is 
more pronounced in signature when looking at the CO SLED of 
M 83, the line ratios of which are above our model galaxies for all 
Jup > 2. Our models MW1, MW2 and MW3 only agree in CO 
luminosity with M51 at the J up = 4 and 5 transitions. 

As for brightness temperature ratios (see middle panel) of 
our simulated MW-like galaxies only agree with the MW within 
the observational errors at CO(4—3), and other wise display a CO 
SLED shape noticably different from the MW. Where our model 
galaxies peak at the CO(3—2) transition in velocity-integrated 
intensity (see bottom panel), similar to M51, the MW and M83 
peak at CO(4—3). However, nowhere do the simulated line ratios 
differ by more than a factor of 2. 


Table Dl. Properties of the three z = 0 MW-like galaxies (MW1, 
MW2, and MW3) used to benchmark SIGAME. 



SFR 

[M 0 yr- 1 ] 

M, 

[10 10 M 0 ] 

Msph 
[10 10 M 0 ] 

/SPH 

Z' 

Rcut 

[kpc] 

MWl 

2.2 

0.97 

1.17 

55% 

2.07 

20 

MW2 

4.1 

3.12 

1.27 

29% 

3.31 

20 

MW3 

6.2 

3.83 

1.69 

31% 

3.90 

20 


Notes. All quantities are determined within a radius R C ut, 
which is the radius where the cumulative radial stellar mass 
function of each galaxy becomes flat. The gas mass (Msph) is 
the total SPH gas mass within R C ut- The metallicity 
( Z' = Z/Zq) is the mean of all SPH gas particles within R c u t- 

The molecular gas masses of MW1, MW2 and MW3 are 
4.6, 5.4 and 8.5 x 10 9 M©, respectively, significantly above (1.0 zb 
0.3) X 10 9 M© for the MW (Heyer & Dame 2015). The corre¬ 
sponding global (R < 20kpc) aco factors are 13.6, 7.7 and 
7.0 M© pc -2 (Kkms -1 ) -1 . For the inner disk (R < lOkpc) of the 
MW, aco,MW — 4.3 zb 0.1 M© pc -2 (Kkms -1 ) -1 (e.g., Strong 
& Mattox 1996; Dame et al. 2001; Pineda et al. 2008; Bolatto 
et al. 2013). However, typically the measurements of qco i n 
the MW do not probe the central (< 1 kpc) region. Excluding 
the central R < 1 kpc region and going out to R = 10 kpc, 
the aco factors of MW1, MW2 and MW3 are 11.6, 6.9 and 
6.2M© pc -2 (Kkms -1 ) -1 , i.e. factors of 2.7, 1.6 and 1.4 above 
aco,MW, respectively. 
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Figure Dl. Global CO SLEDs of our three model galaxies MW1, MW2 and MW3 shown as red (dashed), green (dash-dot) and blue 
(dash-dot-dot-dot) curves, respectively, selected to represent the MW in terms of stellar mass and SFR. The SLEDs are given as absolute 
line luminosities in units of Kkms -1 pc -2 (top panel), as brightness temperature ratios normalised to the CO(l—0) transition (middle 
panel), and as velocity-integrated intensity ratios normalised to CO(1—0) (bottom panel). The observed global CO SLED of the MW 
is shown in grey (Fixsen et al. 1999), and global CO measurements of the local star-forming galaxy M51 are displayed with triangles. 
The CO observations for M51 are from Vlahakis et al. (2013) and Hughes et al. (2013). Also shown, with long-dashed orange lines, are 
observed CO SLEDs of local galaxies with IR luminosities within 20% of that of the MW (Kamenetzky et al. 2015). 
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